{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" }, "colab": { "name": "604_Boundary Value Problem Example 2.ipynb", "provenance": [], "include_colab_link": true } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "view-in-github", "colab_type": "text" }, "source": [ "\"Open" ] }, { "cell_type": "markdown", "metadata": { "id": "YcFC-LNyzHCa" }, "source": [ "# Finite Difference Method\n", "#### John S Butler john.s.butler@tudublin.ie [Course Notes](https://johnsbutler.netlify.com/files/Teaching/Numerical_Analysis_for_Differential_Equations.pdf) [Github](https://github.com/john-s-butler-dit/Numerical-Analysis-Python)\n", "\n", "## Overview\n", "This notebook illustrates the finite different method for a linear Boundary Value Problem.\n", "### Example 2 Boundary Value Problem\n", "To further illustrate the method we will apply the finite difference method to the this boundary value problem\n", "\\begin{equation} \\frac{d^2 y}{dx^2} + 2x\\frac{dy}{dx}+y=3x^2,\\end{equation}\n", "with the boundary conditions\n", "\\begin{equation} y(0)=1, y(1)=2. \\end{equation}\n" ] }, { "cell_type": "code", "metadata": { "id": "Z-AVtF4GzHCe" }, "source": [ "import numpy as np\n", "import math\n", "import matplotlib.pyplot as plt\n" ], "execution_count": 1, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "GRgxWLgVzHCg" }, "source": [ "## Discrete Axis\n", "The stepsize is defined as\n", "\\begin{equation}h=\\frac{b-a}{N}\\end{equation}\n", "here it is \n", "\\begin{equation}h=\\frac{1-0}{10}\\end{equation}\n", "giving \n", "\\begin{equation}x_i=0+0.1 i\\end{equation}\n", "for $i=0,1,...10.$" ] }, { "cell_type": "code", "metadata": { "id": "Fmm-6dldzHCh", "outputId": "a8b8c678-0ac6-4b0c-9c02-43413d3c1473", "colab": { "base_uri": "https://localhost:8080/", "height": 315 } }, "source": [ "## BVP\n", "N=10\n", "h=1/N\n", "x=np.linspace(0,1,N+1)\n", "fig = plt.figure(figsize=(10,4))\n", "plt.plot(x,0*x,'o:',color='red')\n", "plt.xlim((0,1))\n", "plt.xlabel('x',fontsize=16)\n", "plt.title('Illustration of discrete time points for h=%s'%(h),fontsize=32)\n", "plt.show()" ], "execution_count": 2, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAswAAAEqCAYAAAABJAYYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZwlZXno8d8DIyiLgICIbKOCC7ghHRE1FxQE1ABexChXI+NGjDHXBUzwGnVEibuoCYmiGBCIiOAyiEgQHTUo4iAqouwOmyiDLIKIiDz3j7cOXdOc857T3af7dDe/7+dzPn2qzttVT+1PvfVWVWQmkiRJkrpbY9QBSJIkSXOZCbMkSZJUYcIsSZIkVZgwS5IkSRUmzJIkSVKFCbMkSZJUMdSEOSKOjYhsPrv1KLO8VWbxMMev6YmIxa1ls3zU8cxHEfH4iDg6Ii6OiN+35ueMP7+xNa6VPX7frVXm2JmOR3NXRKycrfVSw+H+ubeI2CQi3hMRKyLiloj4c79cZC4bJJfS7Fs0aMHmILxN0/mszFw+EwHNVRHxRmBDgMxcOtpoBhMRLwCe3HQem5krRxjOghcR+wCnAGuNOhZNX0QsARY3nR/NzFtGF01vEfFk4AVN5/L7275ZMysiljZfb8nMj44ylm4i4hHAOcDmo45Fq4uINYC/Bg6k5CKbAb8DfgksAz6dmb8Z8jg3AXZqfcaArVtFppy/DpwwizcyfsKwdIRxTMYLgIOa78uBlSOLZIGLiAcAn2Q8Wf4+8G3gtyMLStO1BNi1+X4sMCcTZsqB6J2t7uUjikMLU2fdugqYcwkz8AHGk+WVwJeA3wB/bvpdMYKY7vciYkvgJOAZE37atPk8FXhzRLw6M780pHG+FviPYQyrGxNmaTj+gvGd9vcz8+mjDEaSFrqIWAQ8r+m8DRjLTCspRiwiNgL+G3hc0+sPwBeBS4CNgP2ARwIPAU6OiH0y8+tDGPUDu/S7ren/gOkO3IRZ92qabMSo45in2pd8vjuyKCqay1AuX5GZi0cdgybH/XNXmwDrNN9/YbI8Z3yI8WT5YuC57SahEfEW4OPA6yh56GcjYtvM/N00x3s75cru+a3PpZQmINtU/m8gJszScLTPbP8wsigk6f7D/e4cExGPpTRnA/gTcMDE+6cy888R8Q/AkyhNNjYFDgXeMZ1xZ+angU93iWk6g11tBAN9KG2Dsvns1qPMsQOUWd4qs7jL74tbvy8fIK4lrfJLK+U2Ag4BzgauB/4I3EVpY/ojSvvTA4ANJ/xfDvhZOeH/lrZ+W9L02wZ4N3ABsAq4Z+I0AmsCuwPvb+ZVJ9Y7gKuBrwCvAtaqTGt7Hvf7LJ7GvH8o8HbKDRe/aebnDZT2u0uBh01mvWr1eyZwAnAlcGezjL4N/C2waNB1dhLr9iOA9wIrgBub+X098C3gLcAGPf6vPb/6fbpuD5OI8QDg9CauOyntCb8EPL/LurqyxzB2a5U5tjKuNYGXUi6h/RL4PXA3cCvwc+Bk4DXANgPEvRHwZuAM4BrKQe3OZl3+WjN/txpkP0C5pPYK4OvN/981cR1u/f+6wOubcVzdjPdW4BfAJ4CdBlkn+30qwwjghc16fDnlRpc/NMvt88D+QExznVg6ifVvSb/tboD5vwh4NfDNZj38A6X26Kgey/CJwKeaef574CbK/veASU7n1s20ngP8ulnuNwLnAe8BHj7EfcFq2xDlSVIvBc4ErqPsG64FvgDsOclhPxM4uplnt7bWhy8CLwPW7PP/91kmXcosaZVZ2vR7MOXYdx5lX3onZd/6KeAxA4yr3+c+sQAPouyvv9bMrz9QkqebgZ8CnwX+Bthsistp5YCxLe3x/1Pa5/dZV4Kyn/4SZb95J1Pc99Mll6I0YfgQZR98O6WpwU8oOcVGw9oGhrANvbsV+3/1KbtXq+wvZzCmlRPn55SGM8wRdlvIXcosb5VZ3OX39oa6fIC4lgywcexKSVAH2cA+NOF/B91prJzwf0tbvy0BXkw5aFZ3NpSkYpDxXQ7sMMA87vdZPJV5D7yyx/S0P7cDB09ivQrgg5QTiZ47Z2C9IW5Ib6PsLGvTsQrYt8+62u8zpY2Ucrnxq32GfRwlmey6LraGtVurzLE9ymxBOaEbZJpW9In9dZTEoN9w/txvP0BJmn7Ybx1u/nd/yglcbZz3UBLnB/RZJ6ufHtO9PeVA1u//z2GKCUOXfUy/z5Je09hj2BPn/2ZNvL2Gfyuwc+v/D6e+HR81wPQF8C7GE49enzvos5+ZxDy9dxsCNqAkyrVx/yf9E911KSeZ/ZbRhcCjK8NZbZn0KLOkVWYppQbv8so4/wi8cJr7tuUT/ncHSsI4yP+eMsXltHLA4S/t8r9T3udX1pUNKSfy3Ya12xSm79j2/wMvoSTIveK9FnjsMLaBIWxDP23F9YI+ZdegnER1yu84QzG115dJL4/OZ8E3yWju1DwNWL/pdS2lMfpKyhnvBsCjgV3o/liatzR//x+ltqzdr+3WShjPoOzIFlHa1HyLspJsSTn7b9u0+XsF5ez3CsYbrT8W2LuJ+VHA2RHxpLzvY1n+g5JkvZjySBUoyUG3u4VvqsTdVUS8HvjXVq9rKI+I+RXwMGAfyg53XeCTEbFODvY4osMpl2XuoNSmXtT034VyJgrl5OcjwMGTjXuiiPgQpeal4xJKjciNlKsB+1EShU2AL0bEizPz1Fb5mxhfF8Yo8xvgLMo61jbpO7WbR/J8GXhOq/d3ge9QatmeQJnXL6fs4KelNb7OowhvoxwELqbUEK5LWa5jjLdP6zWsjwBvavW6gZJ8XEnZ7jan3Cj5F/R/HvwDm7h2bIbzVcoBeT3KtpWt8b4W+HfG23r+jHJ14leUJ5jsCDyXcoLxt5SbTv56wviOoGxjf0ep1QH4F8o2WxURT6Us+w2aXtc33VdSEshtKcvsIcDTgf+JiLHMrO0/evlvyklpv3UPysnGVD0AOJUS79WU+X8dZds4AHg4ZT92ekQ8EvgHypWnOynb8YWUZfws4C+bYb4uIr6XmSd2G2GUa6gnAP+n6XUPZb3/IWW727AZ1tMptZmfjIhFmfnv05jOiY4B9qRsB1+mtIVcj7I9PqUps4Syri3pMR1rUZbJLq3e3wb+h7IN70BZHx4EPB44JyJ2yczLhxD/lpTt92GU7eAsyn5iC8pJ5eaUbeL4iPhxZrb3Ue192webvzdTtoOJrul8iYj1KFeTtmp63djEcAWllnl9yrFrZ0ot71R1ttGNKMdmKNvYxCckfK/dMYR9fjcBnEg5Rt1KWecvoczbMcaf1jFVewL/SNmGvgmcSzlGPgZ4EWX/uAVwSkTsmJl/mub4piwi1qZUGHT8T618Zt4TEd9j/MbNHSkVNnPTMDN05mANMyUJ6/x+ND0u6VNW+l2Av+o3/QPOr6Wsfgb4Bwa4FAl8GNi+8vv6lMtZneEeUynbd3lMdt5TkrT22fm/MqF5COUA+9FWmbvocebIfWsKvkeXS6yUA3OnxurPwBZTPUtshrfXhPEeBqwxocx6lMvnnTI39xpvv/VwijG+vjXMO+lytk7ZaXaSsXtrPHoMb7dWmWO7/P7s1u8/BDauxPZI4NU9fnt5azj3UGoJH9ij7MOB9/RZF++NGVi3EtPTGG+mcQPwvB7lHsXqtSCv6FFueavM4gGW10aUS+xJacJyKF32N5RtuL1eHTfN9WTS6x6Tq2HufI7kvtv6+qxe839SM+0/6TbPKLV7nbIXV+I7pFVuBT1qzyjJxC1NuT/So4nBJOZle71NSnLSbX/0asp+qFO+a00a5bJ/p8xtwF49tqWftcr9gC7NdZh8DXNSTk5f02MdbF8x+OQA82TlAPPvla3yXwXWqZR9InDgNJdX33nSKjvsff7E7ePrwCbTmZ7WsI+dMOxfA7t0KfdYyr6uU67n/GQSV84G/HTbvp/UnncDTuvHWv/z4WHMvz7TvtuUhzPMETI3E+bTWr/3TAAmM/0Dll86YeXqmlxMMZY1KO2uk5KId22iMMjymOy8Bz7XKrOsz/C+2Cp76gDr1TVMaEM+oexJrbJ/P815+P3WsD5eKbeI0vavukH3Ww+nEN8iSi1eZ5h/Wym7PeOJYs8DG/0T5naScp/LtAPGvfaEuN86xeG018Wk1C5W2/xSat+TkiSM9Sm7BeOXOC/rUWZ5a/yLB4i53XbvjQMs306ieTcDtAevDGvS6x6TT5i/WBnWMyaUvYXeScYalJraTtn7VA5Qao87zb1W0qc9KaVWsDO8T091PjbDak/HKioJEKvv43/U5feNKFdmOmX2rwxrS1Zv4rZPn2WyfIB1IYFDK+N8XHtaB5gnKweYf//aKt/zPoFhfQaZJ62yw97nt+fzZcCDhjhdx7aG/WdazZ26lH1tq+wXKuVWToh5up/FXcbRPim5cMBpPbT1PyfO0HrSnvbdpjqcBd8kg3IDU8c6jOZFEtcAnxnWwLJcxjiZcvnigZTL2t8a1vB7iYh1KTcydfxTn3/5J+B/N9/3i4gNs/62tA/3+f0Uxi8979hn3D1FxKMptZFQam7f2atsZt4dEf+PcjkT4OURcWg2W+EMejal9hVKreWnKjH+PCKOp9TuTMfEbWUqXsB43JdRbl4dhsNr8zwinkS5qQrg5MxcURtYZl4XEZ+j3Ly4bURsn5k/n2pwTXOWv2s6r6U8Mqk2/rsj4kjKpdw1gedTmpLMVbVt5JyIuJnxJmufzszrepS9JyK+Bryh6bUj5SamtoMYb0J3RPZprpKZX4mIy4DtKM0bhuXIzLyx8vsHKM2OHgzsGBFPyMwLW7+/iPHt6LzM/GKvAWXmtRHxcUoNPJR5cNrUQwdKwt9zPczMX0TERZRmIZtExJaZee00xzmMfcjQzcI+//2ZOVNP6TgtM39Q+f0Uxpui1I6LnWYsw9KtOef6re+/H3A4d/T4/znn/pAw/4TSZhHgmIhYkpm/muUYzsjMeyb7T80rP59AaVO1HqvvjNobxmOYhYSZ8maezsO/L8rMX9QKZ+ZlEXEBJdY1KTus2sPJz+gz/ktb3x/ap2xN+81DZ2dmv7apZ1NOtDamtG17DKVd70xqx/iVAdafU5l+wvyT1vd3R8TPM/P8SQ5j99b3z05lve/iD/Rfv9vjPXPA4f649f0vuG/iNhlPpqwfAGcNON0Txz9XXT8hEezml4wnzN3aULe128o+rMvvU12W2wEPjYhtMvOqAf+vpvr2scy8IyK+zngb+KdT2mt3tLfhLwwwvpMZT5gnvh1tKs7OzLv6lLmUkjBD2adON2Fu70M+1rQBvmyawxyGmd7nnz698Kqqx8XMvLF1wtrzuJiZPStdhuhBre/91r2OO1vf58xJVjf3h4T5E5S79R9MuVnjqoj4LuUAfC5wbmbeNsMxDJxcRcSalBva/i+lfdIgNpxKUFOwXev7oA3zf8R4cr8d9YS530GuvZzWG3D83UxqOjIzI+LHjB/It2PmE+ZHt77/uGepcT/pX6SvsyjzY0fKDTArIuKnTf/vUbaVfieb7Rs+zhtCTFCaTPS7ceZJre/HRcRxkxzHpv2LDDz+V0TEK2Z5/DPp6gHK3D6J8u2y63b5vT0vr5rCM1Q3pf++pJ87KTdu9fNTxhPmx0z4bbL7y4soScZawMMiYv1pHpsGmQfD2qd2/Bfwz5QmJjsCl0TEeZQE9PuUt6CO4irvTO7zf5eZ108zvppBl+NGDGcZTke7ln2tAf+n/SztO3qWmgP63Z0+7zU1DXsxXquxiHK39uGUmpCbI+K7EfGa5o7mmTDQHfDNHaanUS7NDposQ/fXQc6EjVrfa5cq29rlHlIrmJl31n6ntD/qmM66O6PTMSTtGAc5wAw6HT01taLPA77R6v1EStvmU4HrIuIXEbE0InoleBu3vt8w3Zgag2w/G/cvUtUtcZuMUY9/JvXbLmH1bXO62/FcmJc3D3iVoL3dbTTht0ntZ5qTwnbN53T3M5NdbtPOB5oEfw9KRQmUm+l3pjzJ4jRgVUScHxGHRMRsXn6fyX3+VJ5wMxmTWY6jfhNk+wRs0O2wXas805WX03J/qGEmM89t3j6zD7Av8L8Yf1zUmpS2j88E3hIR+/VrajAFg16WfjvjzUdup7xM5UzKZbNVwJ2dnXhTgzW0dtESQGb+GnhOROxCqTnblZI0d5oDPZbS/u/NEXFQZlYvWw/JINtPe1/2n0y+ecX3+hcZePzfYfLtT6/pX+R+oz0v/5ny9IvJuHKIsWiSMvOSiBijJM77U463j6Mkc0F5JN9TgEMj4kWZWX302DwwjGZnMy4iXsNw2zAfnfd9lfWvW9+3HHA47XITH5E7p8zFhLl9xjvI2dJAbV4y825Km7QvAUTEwykb8j6UG9nWplx2+WpzA9Bkd9LTEhGLGL9p6G7gWX1uXBrmij+odu3HoLVAm7S+T/qZzzNkPkzHZGPcpH+RwWXm9ymXUImIB1PaZ+5NeS7upjSPRmue+3lR61/bteHTaWc+We3xnp09nu87S+O/NDM/NMvjX0h+y/iNo8f2uoFwhm0UEWsMUMvc3u4mtoud1DbcNMdr14TOlf3lpDU3yJ3VfIiIjSnPzX4e5cbtB1Par58WEY/OzGk/R76P+bDPn2lvozS1G5ZTKE92abuEcgKxBrBhRGzS58ZZKM+m77ioZ6k5YC42yWi3bxukPc7WUxlJZv4qM0/KzJdS2sx1NohHMl7LO5sezfhln+/0u8uf8pD72da+eWPQp1S0y13as9TsmtR0NC9RaLernI3paI/jST1LTa7MlGTm7zLz65n5Rsr20blj+wGUp0y0tXd4T52pmLpoXxWasXkxh8e/kMyFeflAVr+PoJcntr5PbPM82f3l9oy3+/z1LNxbM2sy87eZ+eXMPJgyXztXATakvHp8ps2Hff681zSrbF/dq9682jxd6OmtXnP3pSXMzYT5FsYvwW3bzNCavfr83ldmXkK5Q7lj4s0bUJ7tCtxbEzBsA9csNG2d/2qAYbbf+DOMmM9rDfPxEdFtPt0rIh7F+M7pz4wnWqN2Tuv7HhHRr7b+WYzXNtzI7Ow8200E9ov+dz7tP5PBdGTm7ZQXWHRMXAfObn1/+QDb77Cc1fr+ouaKzXRNZvv5AePt73ZqHmM1W4a9nY9ae1keOLIoxh+J2VVEPIhy1aVjYrOe9n7mgAHG96Ie/ztqdzd/h7JuZXkz7dGtXtXjyJDMh33+jMrMxZkZQ/ys7DGqr7S+v6hHmY7dGc99rsrMH9UKj9qcS5ibGx9+2nQ+mNUfMbSaiNif8df4DlO35ym2G/bPxE1f7Uu6T+6TaLyV8qi5foYac2b+nnLzV8d7+/zL+xhvVvPlPs9YnjWZeSnlCSlQapKW9irbnBy1Xwd73Cw8gxlK4tm583ox5c1iXTXt8/9mFmLqZuK28mXKa6ihNHH6x9kIonlOaef1z4sZfzxXX5WTkYG3nyyvo/1E07kGcNRkTqwHOCGqmel902z7DOPPcD0wIp496D9Ocz5O9KaIqM3Pt1COUQA/7vLovS8wftf/0yJiv14DapoIvqHV69hJxjqTOuvXxJsah2Wmnl98r3myz18oTmS8XfeLIuJx3Qo1Oc7bW72On+nApmvOJcyN9gPePx4R90kOI+L5DLBTiYhvRsSbIqLb8z47ZZ4JvKzV6ztdirUvtz2r33in4FLGnyqwLXDExKQ5ItaMiH8C3sHqbb17mYmY/4Xx5yv+74j46MSni0TEAyLiw4zXqvwJeM+Qxj8sS1vf3xgRb5l4sG1e1HI85S5vKFc/2rWrM6Zpc/++Vq+PR8R9XsrQ1GSexvjzsaesWZYfiogdKmW2oNyI1bHattI89/WtrV7/0jxVo+uTXCJii4g4Yjpxt7yJ8drWpRHx4ab9dVcRsXVEvJXez/qd7PbzfsZv3tsDOCMiFlfGv0FEHNQ8q3w6T4Zox7nrLNbqz4imPeu7ms41gWUR8craCUhEPCUijmJ4L8pJSlv9r3Y7djQ3Xb+j1evd9xlAedbvx1q9PhsRe3QZ1mLKs3Y76+p5zOxzfSers36tGxHVZlYRcVJEvCPKOwR6ldme8tjUjm7H25mwtPV9zu3zF4rmoQmfbTrXAk6NiNXaTjfb8sco7dqh1OL3vO8jIpZHRDafpcOPejBz8aY/KDU1b6DcFPBY4BcR8QWa16RS7tx/GqUW4kjKgbKXRwIfAT7UHJguoLy69y7KTUk7M/4GIIDPZ2a3596eTrnhCeDTEbErcDnjB+jfZeZnu/zfQJo3YL0f+HDT6zBg34g4m3Ln6eaUt4E9gjLdn6A87qvm64w3wH9tRGwCrGD1Zx0eN5m2cpl5YUQcQnkFKpTl9IKIOI1SI7oZ5UbK9g7zH3vM05HJzDObpL4zDz9AeX7uGZTa/q0pr9ztHCz/THm9+WzegPRvTQzPptSKLIuI71AOMHdRXmqzL+WG1fa0TNWGlDeMHRIRV1AO3L+kNDV4CGVb3Jvx5PwK4NMTB5KZn42InSgHxaA8VePvorzk4UrKJd7NgTFKO+dgEjXCvWR549zfUt6KuCbwZuDVzTb082Y61ge2otyl37kPoNeNJqdTHocF5UUujwN+RuupDZn5b63vv21qEb9BmV/PAS5vltkKSlOrB1K2kSdSpn/aj7JsXhDUedPdDsB3m/W4fUXnG5k5088OH5rM/GBz5eSVlMdTHQO8KyK+QTkO3EU5FmxHeenLVs2/HjWkEK4Gzqc0dbokIr5MqdRYj3IyNNYqe0L2fovfUsor6XehJMRnRcRyShOBP1KW176Mv/DhRuClc6xG83TG25meFhEnUp4N3Hk2+nWtp+U8jHJT37uivEVwBWVe/oHSxOFJlJPPzknd95j+Gw0HMk/2+QvFoZR1/jGUJ6T8PCK+SDn52ogynx/VlL0bWJJ93ug5qIg4pUvv9mNQ3xURE28yPTkzT6afHOK7uFn9/ee9yixvlbnPu8hb5Xah3Nna6z3mqyjNNZa0+i3tMpxLK8OY+DkOWLtHPIuAb1f+d+WE8ktbvy0ZcB4H5WBfi/EGyoG4Ot2tYR7RZ3iLW2UXt/ov7xPrqygJSG3YtwMHD7peDTB/Bo5vEuv12ykHrtp03Ajs22c4Ay2PKcS3HqX2qRbf8ZQktuu62BrWbq0yx3b5/ehJbCvnAVv3if1NzTrQb1h/GuayppxgXDmJaTmpMqwTa/9bWU9r+4qJn0uB9ae5njyPkkT2GseSyWx3k53/DLhfn+y2Qjkhv3XA+Xgn8H+nOR/v3YYoCflZfcb5WWBRn2GuS2me0S/+nwGPrgyn7zKZzLxtyh/bKr9bjzLrU27E7BX38lbZ/57Een8GsNE0l9ek1tPmf4a1z793XZnONEx1uUwov7JTftixTGMatqacENXm803ACwcY1vJB1+tJrH856DA7n7law0xmfr+5dHMo5WCwNeWM7ypKo/KjMvP6iFjSZ1BPpiTWuwE7Uc5qNqEkwLdRDqzfo7zGt+eTKbK8X/45wGspN4PsQKmRm/al8NY4EnhNRHyF8ra/nSlnY7cwPt2fyszfDDDdnWG+LcqblpZQpn9ThvCik8w8pqlVfi2l1nFbyvy4lVLz+HXgE1me6ztnZea7mxqTg4E9KY/dWZ+yIV9MqV05Ood09juF+G4HnhsRLwJeQVmGG1JOnH4EHJOZywCG0HzztZSTxmdTrro8hlIbug7lqsS1lFq3LwDLmvW1FvuREXECpf31Xs3wNqZsx7+hJAjfAE6abuATxvvNpqnKAZSrMjtTriatR7k6cx2lxvk7wFczs/bc3pdRkoADKbXCG9OnVjjLzTC7RsT/otz08kzKs0Y3pByoV1FqWr4PnJGZ034jYmZ+rblc/g+U2sCtKMtt1C8ymJbM/FiUtzYeRKnZfSJl//0Ayv57JeWel7OB07P/K48nM+5bI2Ivyv0BL6Ps8zemJFPnUvYLfV/dneXejxdFxF8CL6c8znRzynp0A2Wb+iLwX9n/jZazLjNvi4idKevW8ynb8YPpfoX6+ZTpezbl6tG2lH3IWpST56spN8h+LjO/OfPR39dc3+cvFJl5ddPc9cWU/eeTKevCbZQrl8so+cyczhHaos8xT5Kk+4WI6BwQr8rMxaOMRdLcMq9vDpEkSZJmmgmzJEmSVGHCLEmSJFWYMEuSJEkVJsySJElSxZx9rNx8sMkmm+TixYtHHYYkaYjWWmutbcbGxnyElBac888//8bM3LR/SU1kwjwNixcvZsWKno9uliRJmjMi4qpRxzBf2SRDkiRJqjBhliRJkipMmCVJkqQKE2ZJkiSpwoRZkiRJqjBhliRJkipMmCVJkqQKE2ZJkiSpwoRZkiRJqjBhliRJkipMmCVJkqQKE2ZJkiSpwoRZkiRJqjBhliRJkipMmCVJkqQKE2ZJkiSpwoRZkiRJqjBhliRJkipMmCVJkqQKE2ZJkiSpwoRZkiRJqjBhliRJkipMmCVJkqQKE2ZJkiSpwoRZkiRJqlhwCXNE7B0Rl0TE5RFxWJff146Izze//yAiFk/4feuIuD0iDp2tmCVJkjR3LaiEOSLWBI4CngtsDxwYEdtPKPYq4ObM3BY4Enj/hN8/Apwx07FKkiRpflhQCTPwVODyzLwyM+8CTgL2m1BmP+C45vspwO4REQAR8QLgl8BFsxSvJEmS5riFljBvAVzT6r626de1TGbeDdwKbBwR6wH/BLyrNoKIODgiVkTEilWrVg0tcEmSJM1NCy1hno6lwJGZeXutUGYenZljmTm26aabzk5kkiRJGplFow5gyK4Dtmp1b9n061bm2ohYBGwA/BbYGTggIj4AbAjcExF3Zua/zXzYkiRJmqsWWsL8Q2C7iHgEJTF+CfB/JpRZBhwEfB84APhmZibwl50CEbEUuN1kWZIkSQsqYc7MuyPi9cCZwJrAZzLzoog4HFiRmcuAY4DjI+Jy4CZKUi1JkiR1FaVyVVMxNjaWK1asGHUYkiRJfUXE+Zk5Nuo45iNv+pMkSZIqTJglSZKkChNmSZIkqcKEWZIkSaowYZYkSZIqTJglSZKkChNmSZIkqcKEWZIkSaowYZYkSZIqTJglSZKkChNmSZIkqcKEWZIkSaowYZYkSZIqTOpZYkUAABHmSURBVJglSZKkChNmSZIkqcKEWZIkSaowYZYkSZIqTJglSZKkChNmSZIkqcKEWZIkSaowYZYkSZIqTJglSZKkChNmSZIkqcKEWZIkSaowYZYkSZIqTJglSZKkChNmSZIkqcKEWZIkSaowYZYkSZIqTJglSZKkigWXMEfE3hFxSURcHhGHdfl97Yj4fPP7DyJicdP/ORFxfkRc2Px99mzHLkmSpLlnQSXMEbEmcBTwXGB74MCI2H5CsVcBN2fmtsCRwPub/jcC+2TmE4CDgONnJ2pJkiTNZQsqYQaeClyemVdm5l3AScB+E8rsBxzXfD8F2D0iIjMvyMxfNf0vAh4UEWvPStSSJEmasxZawrwFcE2r+9qmX9cymXk3cCuw8YQyLwR+lJl/nKE4JUmSNE8sGnUAc01E7EBpprFnj98PBg4G2HrrrWcxMkmSJI3CQqthvg7YqtW9ZdOva5mIWARsAPy26d4S+BLw8sy8otsIMvPozBzLzLFNN910yOFLkiRprlloCfMPge0i4hERsRbwEmDZhDLLKDf1ARwAfDMzMyI2BE4HDsvMc2YtYkmSJM1pCyphbtokvx44E/gFcHJmXhQRh0fEvk2xY4CNI+Jy4M1A59Fzrwe2Bd4RET9uPg+d5UmQJEnSHBOZOeoY5q2xsbFcsWLFqMOQJEnqKyLOz8yxUccxHy2oGmZJkiRp2EyYJUmSpAoTZkmSJKnChFmSJEmqMGGWJEmSKkyYJUmSpAoTZkmSJKnChFmSJEmqMGGWJEmSKkyYJUmSpAoTZkmSJKnChFmSJEmqMGGWJEmSKkyYJUmSpAoTZkmSJKnChFmSJEmqMGGWJEmSKkyYJUmSpAoTZkmSJKnChFmSJEmqMGGWJEmSKkyYJUmSpAoTZkmSJKnChFmSJEmqMGGWJEmSKkyYJUmSpAoTZkmSJKnChFmSJEmqMGGWJEmSKkyYJUmSpAoTZkmSJKliwSXMEbF3RFwSEZdHxGFdfl87Ij7f/P6DiFjc+u2tTf9LImKvviM7/3xYvBhOPHGYkzB6J55YpmuNNZy++cjpm78W8rSB0zffOX3zVzNtO8FOow5l3srMBfMB1gSuAB4JrAX8BNh+QpnXAZ9ovr8E+Hzzffum/NrAI5rhrFkb306QCZnrrJN5wgm5IJxwQpmezrQ5ffOL0zd/LeRpy3T65junb/5qTdtOkDkH8rX5+Bl5AEOdGNgFOLPV/VbgrRPKnAns0nxfBNwIxMSy7XK9Pju1N6ytt87cddfM44/PzMz8/e9L90knle5bbindp55auletKt3LlpXu668v3WecUbqvvrp0n3VW6b7iitK9fHnpvvji0n3OOaX7wgtL93nnle4LLijdF1xQus87r3RfeGHpPuec0n3xxaV7+fLS/fCH52o7jM7noQ8tv19/fSm/bFnpXrWqdJ96aum+5ZbSfdJJpfv3vy/dxx9fuu+6q3T/53+W7o6jj87cfffx7qOOytx77/Huj340c599xrs/+MHM/fcf737vezNf/OLx7sMPz3zpS8e73/72zCVLMrfZpvv0rb9+5uteN17+DW8on47XvS7zkEPGu1/zmszDDhvvXrKkjKPjpS8tMXS8+MUlxo799y/T0LHPPmUaO/beu8yDjt13L/OoY9ddyzzMLPO0s+71mr611577694VV5Tus84q3VdfXbrPOGN83atN31xf9zoOO6ysPx2HHFLWr17Ttt5682Pdy6zv92rLbj6se5n1/V5t+ub6utdR2+/1mr4NNhgvP1fXvcz++70ttug+fZttVn6fy+teZn2/11p2JsxT/yy0JhlbANe0uq9t+nUtk5l3A7cCGw/4v0TEwRGxIiJWrPbDNddMLDo/XX999/6rVs1uHDPl6qu797/tttmNY6b0mr4//nF245gpC3n6ek3b7bfPbhwzZSEvO7j/Tt+tt85uHDPlV7/q3v+GG2Y3jpnQa9lpUiIzRx3D0ETEAcDemfnqpvtvgJ0z8/WtMj9rylzbdF8B7AwsBc7NzBOa/scAZ2TmKb3GNxaR92bN22wDK1cOfZpm3eLFcNVV9+3v9M0PTt/8tZCnDZy++c7pm79a0zYGrMiMkcYzTy20GubrgK1a3Vs2/bqWiYhFwAbAbwf83+7WWQeOOGJqEc81RxxRpqfN6Zs/nL75ayFPGzh9853TN391mzZN3qjbhAzzQ2mTfCXlpr3OTX87TCjz96x+09/JzfcdWP2mvysZ5Ka/bbZZGDcFtJ1wQpmuCKdvPnL65q+FPG2ZTt985/TNX8202YZ56p8F1SQDICKeB3yU8sSMz2TmERFxOLAiM5dFxAOB44EdgZuAl2Tmlc3/vg14JXA38MbMPKM2rrGxsVyxYkWtiCRJ0pwQEedn5tio45iPFlzCPJtMmCVJ0nxhwjx1C60NsyRJkjRUJsySJElShQmzJEmSVGHCLEmSJFWYMEuSJEkVJsySJElShQmzJEmSVGHCLEmSJFWYMEuSJEkVJsySJElShQmzJEmSVGHCLEmSJFWYMEuSJEkVJsySJElShQmzJEmSVGHCLEmSJFWYMEuSJEkVJsySJElShQmzJEmSVGHCLEmSJFWYMEuSJEkVJsySJElShQmzJEmSVGHCLEmSJFWYMEuSJEkVJsySJElShQmzJEmSVGHCLEmSJFWYMEuSJEkVJsySJElSxYJJmCPiIRFxVkRc1vzdqEe5g5oyl0XEQU2/dSLi9Ii4OCIuioj3zW70kiRJmqsWTMIMHAacnZnbAWc33auJiIcA7wR2Bp4KvLOVWH8oMx8L7Ag8IyKeOzthS5IkaS5bSAnzfsBxzffjgBd0KbMXcFZm3pSZNwNnAXtn5h2Z+S2AzLwL+BGw5SzELEmSpDluISXMm2Xm9c33XwObdSmzBXBNq/vapt+9ImJDYB9KLbUkSZLu5xaNOoDJiIhvAA/r8tPb2h2ZmRGRUxj+IuBzwMcz88oeZQ4GDgbYeuutJzsKSZIkzTPzKmHOzD16/RYRv4mIzTPz+ojYHLihS7HrgN1a3VsCy1vdRwOXZeZHKzEc3ZRjbGxs0km5JEmS5peF1CRjGXBQ8/0g4CtdypwJ7BkRGzU3++3Z9CMi3gNsALxxFmKVJEnSPLGQEub3Ac+JiMuAPZpuImIsIj4NkJk3Ae8Gfth8Ds/MmyJiS0qzju2BH0XEjyPi1aOYCEmSJM0tkWmrgqkaGxvLFStWjDoMSZKkviLi/MwcG3Uc89FCqmGWJEmShs6EWZIkSaowYZYkSZIqTJglSZKkChNmSZIkqcKEWZIkSaowYZYkSZIqTJglSZKkChNmSZIkqcKEWZIkSaowYZYkSZIqTJglSZKkChNmSZIkqcKEWZIkSaowYZYkSZIqTJglSZKkChNmSZIkqcKEWZIkSaowYZYkSZIqTJglSZKkChNmSZIkqcKEWZIkSaowYZYkSZIqTJglSZKkChNmSZIkqcKEWZIkSaowYZYkSZIqTJglSZKkChNmSZIkqcKEWZIkSaowYZYkSZIqFkzCHBEPiYizIuKy5u9GPcod1JS5LCIO6vL7soj42cxHLEmSpPlgwSTMwGHA2Zm5HXB2072aiHgI8E5gZ+CpwDvbiXVE7A/cPjvhSpIkaT5YSAnzfsBxzffjgBd0KbMXcFZm3pSZNwNnAXsDRMR6wJuB98xCrJIkSZonFlLCvFlmXt98/zWwWZcyWwDXtLqvbfoBvBv4MHBHbSQRcXBErIiIFatWrZpmyJIkSZrrFo06gMmIiG8AD+vy09vaHZmZEZGTGO6TgUdl5psiYnGtbGYeDRwNMDY2NvA4JEmSND/Nq4Q5M/fo9VtE/CYiNs/M6yNic+CGLsWuA3ZrdW8JLAd2AcYiYiVlnjw0IpZn5m5IkiTpfm0hNclYBnSeenEQ8JUuZc4E9oyIjZqb/fYEzszM/8jMh2fmYuCZwKUmy5IkSYKFlTC/D3hORFwG7NF0ExFjEfFpgMy8idJW+YfN5/CmnyRJktRVZNoMd6rGxsZyxYoVow5DkiSpr4g4PzPHRh3HfLSQapglSZKkoTNhliRJkipMmCVJkqQKE2ZJkiSpwoRZkiRJqjBhliRJkipMmCVJkqQKE2ZJkiSpwoRZkiRJqjBhliRJkipMmCVJkqQKE2ZJkiSpwoRZkiRJqjBhliRJkipMmCVJkqQKE2ZJkiSpwoRZkiRJqjBhliRJkipMmCVJkqQKE2ZJkiSpwoRZkiRJqjBhliRJkipMmCVJkqQKE2ZJkiSpIjJz1DHMWxFxG3DJqOPQlG0C3DjqIDRlLr/5y2U3v7n85q/HZOb6ow5iPlo06gDmuUsyc2zUQWhqImKFy2/+cvnNXy67+c3lN39FxIpRxzBf2SRDkiRJqjBhliRJkipMmKfn6FEHoGlx+c1vLr/5y2U3v7n85i+X3RR5058kSZJUYQ2zJEmSVGHCLEmSJFWYMA8gIvaOiEsi4vKIOKzL72tHxOeb338QEYtnP0r1MsDye3NE/DwifhoRZ0fENqOIU/fVb9m1yr0wIjIifNTVHDLI8ouIv262v4si4r9mO0b1NsC+c+uI+FZEXNDsP583ijh1XxHxmYi4ISJ+1uP3iIiPN8v2pxHxlNmOcb4xYe4jItYEjgKeC2wPHBgR208o9irg5szcFjgSeP/sRqleBlx+FwBjmflE4BTgA7MbpboZcNkREesDbwB+MLsRqmaQ5RcR2wFvBZ6RmTsAb5z1QNXVgNvfPwMnZ+aOwEuAf5/dKFVxLLB35ffnAts1n4OB/5iFmOY1E+b+ngpcnplXZuZdwEnAfhPK7Acc13w/Bdg9ImIWY1RvfZdfZn4rM+9oOs8FtpzlGNXdINsewLspJ6l3zmZw6muQ5fca4KjMvBkgM2+Y5RjV2yDLL4EHN983AH41i/GpIjO/A9xUKbIf8NkszgU2jIjNZye6+cmEub8tgGta3dc2/bqWycy7gVuBjWclOvUzyPJrexVwxoxGpEH1XXbNZcStMvP02QxMAxlk23s08OiIOCcizo2IWo2YZtcgy28p8LKIuBb4GvAPsxOahmCyx8b7PV+NLTUi4mXAGLDrqGNRfxGxBvARYMmIQ9HULaJcEt6NcmXnOxHxhMy8ZaRRaVAHAsdm5ocjYhfg+Ih4fGbeM+rApGGzhrm/64CtWt1bNv26lomIRZRLU7+dlejUzyDLj4jYA3gbsG9m/nGWYlNdv2W3PvB4YHlErASeBizzxr85Y5Bt71pgWWb+KTN/CVxKSaA1eoMsv1cBJwNk5veBBwKbzEp0mq6Bjo0aZ8Lc3w+B7SLiERGxFuXGhmUTyiwDDmq+HwB8M30jzFzRd/lFxI7AJynJsm0o547qssvMWzNzk8xcnJmLKe3P983MFaMJVxMMsu/8MqV2mYjYhNJE48rZDFI9DbL8rgZ2B4iIx1ES5lWzGqWmahnw8uZpGU8Dbs3M60cd1Fxmk4w+MvPuiHg9cCawJvCZzLwoIg4HVmTmMuAYyqWoyymN7F8yuojVNuDy+yCwHvCF5l7NqzNz35EFLWDgZac5asDldyawZ0T8HPgz8JbM9OrcHDDg8jsE+FREvIlyA+ASK4vmhoj4HOVkdJOmjfk7gQcAZOYnKG3OnwdcDtwBvGI0kc4fvhpbkiRJqrBJhiRJklRhwixJkiRVmDBLkiRJFSbMkiRJUoUJsyRJklRhwixJkiRVmDBLkiRJFSbMkiRJUoUJsySNWESsGxEXR8R5EfGAVv89I+KeiPj7UcYnSfd3vulPkuaAiNgROBc4MjMPi4jNgJ8AP8jM/UYbnSTdv5kwS9IcERFvAj4E7AUcCjwBeFJm3jjSwCTpfs6EWZLmiIgI4HTg2cBawHMy8+zRRiVJsg2zJM0RWWowjgfWBn5isixJc4MJsyTNERHxMOBjwI+AJ0XEG0YckiQJE2ZJmhOa5hjHAX8E9gA+Crw/Ip440sAkSbZhlqS5ICIOAT4APDszvx0Ra1GemrE2MJaZfxhpgJJ0P2YNsySNWEQ8BfgX4L2Z+W2AzLwLOBBYDHxkdNFJkqxhliRJkiqsYZYkSZIqTJglSZKkChNmSZIkqcKEWZIkSaowYZYkSZIqTJglSZKkChNmSZIkqcKEWZIkSar4/2h7kevQ42ybAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "xVUcR9Y2zHCh" }, "source": [ "## The Difference Equation\n", "To convert the boundary problem into a difference equation we use 1st and 2nd order difference operators.\n", "\n", "The first derivative can be approximated by the difference operators:\n", "\\begin{equation}D^{+}U_{i}=\\frac{U_{i+1}-U_{i}}{h_{i+1}} \\ \\ \\ \\mbox{ Forward,} \\end{equation}\n", "\\begin{equation}D^{-}U_{i}=\\frac{U_{i}-U_{i-1}}{h_i} \\ \\ \\ \\mbox{ Backward,} \\end{equation}\n", "or\n", "\\begin{equation}D^{0}U_{i}=\\frac{U_{i+1}-U_{i-1}}{x_{i+1}-x_{i-1}} \\ \\ \\ \\mbox{ Centered.} \\end{equation}\n", "The second derivative can be approxiamed \n", "\\begin{equation}\\delta_x^{2}U_{i}=\\frac{2}{x_{i+1}-x_{i-1}}\\left(\\frac{U_{i+1}-U_{i}}{x_{i+1}-x_{i}}-\\frac{U_{i}-U_{i-1}}{x_{i}-x_{i-1}}\\right) \\ \\ \\ \\mbox{ Centered in $x$ direction} \\end{equation}" ] }, { "cell_type": "markdown", "metadata": { "id": "xrfsYLYzzHCi" }, "source": [ "Given the differential equation\n", "\\begin{equation} \\frac{d^2 y}{dx^2} + 2x\\frac{dy}{dx}+y=3x^2,\\end{equation}\n", "the difference equation is,\n", "\\begin{equation}\\frac{1}{h^2}\\left(y_{i-1}-2y_i+y_{i+1}\\right)+2x_i\\frac{y_{i+1}-y_{i-1}}{2h}+y_i=3x^2_i \\ \\ \\ i=1,..,N-1. \\end{equation}\n", "\n", "Rearranging the equation we have the system of N-1 equations\n", "\\begin{equation}i=1: (\\frac{1}{0.1^2}-\\frac{2x_1}{0.2})\\color{green}{y_{0}} -\\left(\\frac{2}{0.1^2}-1\\right)y_1 +(\\frac{1}{0.1^2}+\\frac{2x_1}{0.2}) y_{2}=3x_1^2\\end{equation}\n", "\\begin{equation}i=2: (\\frac{1}{0.1^2}-\\frac{2x_2}{0.2})y_{1} -\\left(\\frac{2}{0.1^2}-1\\right)y_2 +(\\frac{1}{0.1^2}+\\frac{2x_2}{0.2}) y_{3}=3x_2^2\\end{equation}\n", "\\begin{equation} ...\\end{equation}\n", "\\begin{equation}i=8: (\\frac{1}{0.1^2}-\\frac{2x_8}{0.2})y_{7} -\\left(\\frac{2}{0.1^2}-1\\right)y_8 +(\\frac{1}{0.1^2}+\\frac{2x_8}{0.2})y_{9}=3x_8^2\\end{equation}\n", "\\begin{equation}i=9: (\\frac{1}{0.1^2}-\\frac{2x_9}{0.2})y_{8} -\\left(\\frac{2}{0.1^2}-1\\right)y_9 +(\\frac{1}{0.1^2}+\\frac{2x_9}{0.2}) \\color{green}{y_{10}}=3x_9^2\\end{equation}\n", "where the green terms are the known boundary conditions.\n", "\n", "Rearranging the equation we have the system of 9 equations\n", "\\begin{equation}i=1: -\\left(\\frac{2}{0.1^2}-1\\right)y_1 +(\\frac{1}{0.1^2}+\\frac{2x_1}{0.2})y_{2}=-(\\frac{1}{0.1^2}-\\frac{2x_1}{0.2})\\color{green}{y_{0}}+3x_1^2\\end{equation}\n", "\\begin{equation}i=2: (\\frac{1}{0.1^2}-\\frac{2x_2}{0.2})y_{1} -\\left(\\frac{2}{0.1^2}-1\\right)y_2 +(\\frac{1}{0.1^2}+\\frac{2x_2}{0.2}) y_{3}=3x_2^2\\end{equation}\n", "\\begin{equation} ...\\end{equation}\n", "\\begin{equation}i=8: (\\frac{1}{0.1^2}-\\frac{2x_8}{0.2})y_{7} -\\left(\\frac{2}{0.1^2}-1\\right)y_8 +(\\frac{1}{0.1^2}+\\frac{2x_8}{0.2}) y_{9}=3x_8^2\\end{equation}\n", "\\begin{equation}i=9: (\\frac{1}{0.1^2}-\\frac{2x_9}{0.2})y_{8} -\\left(\\frac{2}{0.1^2}-1\\right)y_9 =-(\\frac{1}{0.1^2}+\\frac{2x_9}{0.2}) \\color{green}{y_{10}}+3x_9^2\\end{equation}\n", "where the green terms are the known boundary conditions.\n", "This is system can be put into matrix form \n", "\\begin{equation} A\\color{red}{\\mathbf{y}}=\\mathbf{b} \\end{equation}\n", "Where A is a $9\\times 9 $ matrix of the form\n", "\n", "which can be represented graphically as:" ] }, { "cell_type": "code", "metadata": { "id": "9LLj2_c6zHCj", "outputId": "1d961608-227a-483b-83f3-9aef987b4c5c", "colab": { "base_uri": "https://localhost:8080/", "height": 297 } }, "source": [ "A=np.zeros((N-1,N-1))\n", "# Diagonal\n", "for i in range (0,N-1):\n", " A[i,i]=-(2/(h*h)-1) # Diagonal\n", "\n", "for i in range (0,N-2): \n", " A[i+1,i]=1/(h*h)-2*(i+1)*h/(2*h) ## Lower Diagonal\n", " A[i,i+1]=1/(h*h)+2*(i+1)*h/(2*h) ## Upper Diagonal\n", " \n", "plt.imshow(A)\n", "plt.xlabel('i',fontsize=16)\n", "plt.ylabel('j',fontsize=16)\n", "plt.yticks(np.arange(N-1), np.arange(1,N-0.9,1))\n", "plt.xticks(np.arange(N-1), np.arange(1,N-0.9,1))\n", "clb=plt.colorbar()\n", "clb.set_label('Matrix value')\n", "plt.title('Matrix A',fontsize=32)\n", "plt.tight_layout()\n", "plt.subplots_adjust()\n", "plt.show()" ], "execution_count": 3, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU4AAAEYCAYAAAAzhB+DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3debhcVZnv8e+PhDECQUAFogIt0iIXGXJDlEE0IJNC4xi8yCDdaa7YDHqvrW23gl5ut31tbbV9RAQElFEQQUEgqIATQwIBAwEJEDAhgAwCgkBy8t4/1qqkqJw6tXfNVef34dnP2bWHtdepkDdr77X2ehURmJlZcWv0ugJmZoPGgdPMrCQHTjOzkhw4zcxKcuA0MyvJgdPMrCQHTus6SXtJiryc1ev6mJXlwNkBkhZVBYaQ9LSk9Uqcf2LN+SHpmE7WueraJ+XlhG5cbxBIOqDmz+KSXtfJesuBszs2AN5f4vijOlWRAj6fFwfOVY6u+fweSZv2pCbWFxw4O6/yalahYCjpvwP/LX9c0ZEaWWE5QL4nf3wh/1wT+EhvamT9wIGz836ef75d0tYFjv9o/rkCuK4jNeqxiLguIpSXI3tdnwYOJwVKgH8G/pLXa1uhNo44cHbemfmngCPHOlDSOsDM/PFa4A+dq5YVVAmQLwCnA5flz9tJmt6bKlmvOXB23q3AHXn9CEljfefvAybn9e8WvYCkjSUdIelsSfMkPSVpWf45X9K3Je06xvlbVjo+qja/fpQOqpB0Xc25Z1Xt2ytv207SV/K1n6ztPW/Uqy7pi1X772rUsSbpa1XHz5W0VqEvrgFJbwXelD9eHhFPA+dUHeJW5zjlwNkdlVbn64C9xziu8hz0KeDSIgVL2g94BDiLdFv5FlLwnZh/vhmYBdyYg9zaZStfhqTjgduAE/O1N2qimJOB3+b1NwFfHeN6BwLH5Y/PAYdGxEtNXHM01YGxEjCvIX3fAB+SNKlN17IBMrHXFRgnvg/8O7AWKTheU3uApNcD78wfz4uIFyUVKbsSJJcDN5Nat4+Sbi03BqYBe5AeFRyRzzmypowngf+d1/9f/vkU8H9Hud5Yjw8+BFSGTV1PCn7PAVsCfyrwuwAQEcslfRiYB2wIzJJ0TUS8bBiQpM14ecv8HyLi90WvMxZJryD9PpC+z6tz3UYknQt8Elgf+CAl7g5sSESElzYvwCJSb3oAf523XZw//wXYaJRzTqo6Z5e87ayqbcfUudaewP8EJo9Rnx2BB6rK2mOMYyvHLCr4u1bXMYAngLc3OGevquPPGuO4D1Ud9yTw2qp9AmZX7b+gzX+GR1eV/dWafTtU7ftVr/9/89L9xbfq3VO5XV8HOLR6h1LT8sj88Y6ImFu00Ii4ISK+FRF1W3QRMQ94b9WmWUXLb8LMiLi+HQVFxIWs+t42As6VNCF//hSrHns8CPx9O65ZZbTb9Eq97gBuzx93k7Rtm69tfc6Bs3uuBh7O6x+t2TcDeH1e78htX0TcBtyTP+7ZiWsAv46I2W0u8zhW1XsP4J/zWNcv5m0jwIcjddy0haQ3AW/NH+fn766WO4nGMT/j7JJIz8bOAT4N7CJph9xygVWB9CXS89Cm5I6f7YFtSc8G1yPd0lZUxiO+TtI6EfEC7XVFm8sjIp6TdChwI+kZ8b8Af8eq3+ULEfGbNl+2bmuzyrnAl0h/hw6X9E8RsbzN9bA+5cDZXd8lBU5InUQnSpoMHJK3/TgiHi9baO4kOYk0BnSDgqdNZlXvcLvc3ebygNRalvSPpN71CcAWedcvgVPaeS1J1W8FrSAFyNHq9Kika4ADgFcD7wZ+1M66WP/yrXoXRerx/XX+eFj+S3oo6bknNHGbLmkn4Hek55ZFgyZV12yntt0uj+JrwIKqzy8Ch0XESJuv8x7gVXn9ZxHx8BjH+nZ9nHKLs/vOBHYDNiH9Ja3cpi8FripTUL41/wFp2BHAncCppGFADwHPVt+OS7qezj3fhM6+W38IqwajA6wN7AOc0ebrVAfAjRtMe1f9j8/+kjZvEGhtSDhwdt9FwNeBSaQOju3y9nOaaD0dCPxVXr+JNAzoxTGO37Bk+X1B0mtJrzvW+pqkX0b7xm5OAfat2rRzXoqYQBon+6/tqIv1N9+qd1lE/JnUSoRVQROa602vflf6m2MFzfxY4I1NXKOn8iuq32fVG0hnsupNoknA+e16xZI0JGxCo4PGUDtawoaUA2dvnFnz+TcRcc+oR46t+nXGJxsceyCwboEyKz3DrQSQdvonVj1e+D1peNKnSW8VQWoRjvaGUyl5LG311H/7xqoZnMZcgHvzOW+Q9PZW62L9z4GzByLil6RXG7+Zl5OaLOqJqvVd6h0kaV2KB5dKB08z75i3VZ5k4/P540uk99Cfi/Qu+qHA83nfJyS9q8XLvQOoTPv3KPCzEueeV7XuTqJxwIGzRyLiUxHx8bw0O2j8hqr1T+WB4S8jaXPgp6SOlajdP4pKy3eSpGlN1qtlkjYgDQWqPIf/bETcWtkfEXezapZ6AWe3OCt7dcC7sOTz5uohS++XNJDPkq04dw4NtqtIt6w7kp733SjpCtLwpJdIsxMdSBoI/3PSAPLdG5R5BfC2vP7jPKHFg6Q3dACWREShmZtadCqwVV6fDfxH7QER8R1J+5Km43sN6b35A8teKI+lrX4lddSxm/VExL2S5gBTSY9DDs31t2HV65flh3FhlEk+miznrKpy6k3ysRWwsOq40ZbZpCFL11Vt27JOeeuTxkvWK+u6Meq4V8Hfa6+qc84aZf8RVfv/CGw2RlkbkWZsqhx/XBPf87FV59/b5J/VCVVl3NLr/we9dHbxrfqAi4gHgJ1InShzgGdJg8MfAn5Mav28KyKeqFvIy8t7FtiVlCbit6ROp669SijpDcB/VW06KiKW1js+Ip4CDmPVGNJ/l7RDyctW36afV/eosV3Aqlb51CbqYANE+V9LMzMryC1OM7OSHDjNzEpy4DQzK8mB08ysJAdOM7OShmIA/FpaO9ah/Vla/2qHP7e9TIB7nt+kI+Wa1Xrp0T+x/JnnC6VLLWrfd0yKJ55s/GLV3DtevDoi9mvntfvFUATOdZjErprR9nIv+emNbS8TYI+5R3akXLNa936i3dOVwhNPjnDz1a9reNyEze4d2haCb9XNrJQgWBbLGy6NSDpT0mOS5ldte6Wk2ZLuzT83ytsl6euSFkq6Q1LReVI7woHTzEoJYAXRcCngLKD2Vv7TpJQl25BmqKrk6Nof2CYvs4BvteN3aZYDp5mVtqLAf41ExA2sPo/swcDZef1s4G+qtp8TyY3A5JyksCeG4hmnmXVPEIx07lXtV1fNTfAIKYMopMymf6g6bnHeVnceg05y4DSzUgJYViwv3yZ5ur2K0yLitMLXiQhJfTmZhgOnmZVW8Bnm4xExtWTRj0raLCKW5lvxx/L2JcBrq46bkrf1RFefcY7Wi1azv696zsxsdQGMRDRcmnQ5aT5W8s/LqrYfnmPEdODpsaYb7LRudw6dxeq9aNX6qufMzFYXBMsKLI1IOp805+u2khZLOhr4N2AfSfcCe+fPAFcC95Mm7f4O8LFO/G5FdfVWPSJukLTlGIes7DkjpYGYXGm2d6WCZtZYwEgbnjxGxKF1dq32NkuOCce2ftX26LfhSPV6zlYjaZakOZLmLKNuOnEza7M0jrPxMswGtnMo986dBrCBXtmXPW9mw0mM0NbX3wdOvwXOvuo5M7PVBbAsxnfg7Ldb9b7qOTOz1QUwkludYy3DrKstztyLthdpYOxi4PPAmgARcSqp5+wAUs/Z88BR3ayfmRWzYpy3OLvdq16vF62yv696zsxsdSsQLzGh19XoqX57xmlmA8AtTjOzEirPOMczB04zK0mMRL/1K3eXA6eZlZJmR/IzTjOzwiLc4hyKwLn1Dn/mgit/0/Zy3zflbW0vE+CXi8/qSLngRHDWHSv8jNPMrLhAvBTjO3SM79/ezEpLk3z4Vt3MrJQRj+M0MysuECNucZqZFZdmRxrfoWN8//ZmVlog36r3ugJmNnjGe+dQt7NcvlbSLyTdJelOScePcowzXZr1sQixLCY0XIZZt1ucy4FPRsStktYH5kqaHRF3VR1TnelyV1Kmy127XE8zqyOlB3aLs2siYmlE3JrXnwUWsHoytpWZLiPiRmByTkxvZn1ihDUaLsOsZ884c5rgnYCbanbVy3TpFBpmfSCQ5+PsxUUlvQK4BDghIp5psoxZwCyAKVsM9/MUs37i4Ug9CJyS1iQFzXMj4oejHFIo02V1euAd37KW0wObdU17krFJWgQ8C4wAyyNiqqRXAhcCWwKLgA9GxFMtX6zNut2rLuAMYEFEfKXOYc50adbHAlgRazRcCnpHROwYEVPz508DP4uIbYCf5c99p9stzt2AjwC/kzQvb/sn4HXgTJdmg6KDqTMOJmXCBTgbuA74x05drFndznL5Kxj7G3emS7P+FiGWrSgUOjaRNKfq82n5EdvKooBrJAXw7bzv1VV3mI8Ar25LpdtsfD/hNbPS0rRyhVqcj1fdgo9m94hYIulVwGxJd7/sOhGRg2rfceA0s5LakzojIpbkn49JuhSYBjwqabOIWJrHbz/W8oU6YLhHqZpZ26XhSK29cilpUn57EEmTgHcB80mdw0fkw44ALuvcb9I8tzjNrJQ2DYB/NXBpGmjDROC8iLhK0i3ARZKOBh4EPtjqhTrBgdPMSmt1dqSIuB94yyjbnwBmtFR4FwxF4Pz985vwzls/2vZyZ//hjLaXCfC+Kbt1pFyAXy8+pyPl7jb38I6Ua4MnwqkzhiJwmln3BGL5ivH9mrMDp5mV1sEB8APBgdPMSkmvXDpwmpmVIJYP+QzvjThwmlkp7hxy4DSzJpSY/WgoOXCaWSmeAd6B08xKCmD5OG9xdnsi43Uk3Szp9pwe+ORRjllb0oU5PfBNOTeRmfWRNk5kPJC6/du9CLwzIt4C7Ajsl2d5r3Y08FREvAH4KvClLtfRzMYS6Va90TLMup0eOCLiz/njmnmpnW/vYNLMzwAXAzNyyg0z6wOVW/VGyzDr+m8naUJOm/EYMDsi6qYHjojlwNPAxqOUM0vSHElzlj/9XKerbWZZZQC8W5xdFBEjEbEjKXvlNEnbN1nOaRExNSKmTtxwUnsraWZjcuDskYj4E/ALYL+aXSvTA0uaCGwIPNHd2plZPZXhSA6cXSJpU0mT8/q6wD7A3TWHVc8A/X7g5zmBm5n1g/Azzm6P49wMOFvSBFLQvigifiLpC8CciLiclHf9e5IWAk8CM7tcRzMbgyf56H564DuAnUbZ/rmq9ReAD3SzXmZWjgOnmVkJgRhZMdy34o04cJpZaQXzqg8tB04zKyXCt+oOnGZWWjhwmpmV4WecDpxj2OfWoztS7tUdSjsMcEiHUg/fuOR7HSl3+pyPdKRc65x2DUeStB/wNWACcHpE/FvLhXbJ+P5nw8zKi/Scs9EyljyW+5vA/sB2wKGStut85dvDLU4zKyWAkdbfDJoGLIyI+wEkXUCaGe2uVgvuBgdOMyupLe+ir5wFLVsM7Npqod3iW3UzK63grfomlakf8zKrx9V+GUm7Szoqr28qaaui57rFaWalFRyO9HhETK2zb+UsaNmUvK0rJH0emApsC3yXNKn694FCvasOnGZWSgTtGI50C7BNbuUtIU3m8+FWCy3hENK8GbcCRMTDktYverIDp5mV1upEjxGxXNLHgatJw5HOjIg721C1ol6KiJAUAJJKzYbek8CZhyLMAZZExLtr9q0NnAPsQprA+EMRsajrlTSzutrx5lBEXAlc2XptmnKRpG8DkyX9HfBR4DtFT+5Vi/N4YAGwwSj7Vma5lDSTlOXyQ92snJnVF+3pVe+piPiypH2AZ0jPOT8XEbOLnt/1wClpCnAgcArwiVEOORg4Ka9fDPyXJHkWeLM+EcPxrnoOlIWDZbVeDEf6T+BTwIo6+53l0qzfRYGlj0l6VtIzeXlB0oikZ4qe39UWp6R3A49FxFxJe7VSVkScBpwGsN42m/f5H5PZcBn0FmdErOxBlyTSne70oud3u8W5G3CQpEXABcA7JX2/5hhnuTTrYwGsWKGGy6CI5EfAvkXP6XbOoc8AnwHILc7/FRGH1RxWyXL5W5zl0qz/BDDgLU5J7636uAZpMPwLRc/vi3GcznJpNliGoCnznqr15cAi0u16IT0LnBFxHXBdXneWS7OBIWKAbsVHExFHtXJ+X7Q4zWzADGiLU9I3GKP2EXFckXIcOM2snMEexzmnHYU4cJpZeQPa4oyIs9tRjgOnmZU3uC1OIM2/CfwjKW3HOpXtEfHOIud7ImMzK2/A3xwCziXNl7EVcDKpV/2Woie7xdkD+3YoeybAlQ+d3pFyD9pi946Ue/OScztS7rQ5/6Mj5RpDMY4T2DgizpB0fERcD1wvyYHTzDon6s00MTiW5Z9LJR0IPAy8sujJDpxmVt7gtzj/j6QNgU8C3yBNcXli0ZMdOM2sNPX/M8xGboqIp0mzr72j7MnuHDKzcop0DPV/YP21pGskHS1po7InO3CaWUmCFQWWPhYRbwT+GXgzMFfSTyTVTjhUlwOnmZU3+C1OIuLmiPgEMI00oVDhwfEOnGZW3oAHTkkbSDpC0k+B3wBLSQG0kF7kHFoEPAuMAMtrE9bn2Zi/BhwAPA8cGRG3drueZlZHgPr8VryA24EfAV+IiN+WPblXverviIjH6+zbH9gmL7sC38o/zaxf9HmLsoCtW5kgvR9v1Q8GzsnT2d9Iynu8Wa8rZWbDo9WsEnUDZ876Ni2vr8ifx1oek3SZpK0b1Rm4RtJcSbNG2b8yy2W2OG+rrZ+zXJr1iKLx0lL50kmSlkial5cDqvZ9RtJCSfdIKpwnqJ3GulX/AiloVdYbfRUbAIeQMk/uPcZxu0fEEkmvAmZLujsibiha4QpnuTTrkaBbw42+GhFfrt4gaTtSOp03A5sD10p6Y0SMtHoxSWtFxEtFjq0bOCPi5Kr1kwpe+Drg/LGOiYgl+edjki4l9WRVB86VWS6zKXmbmfWL3jVVDgYuiIgXgQdybrJppOSOheVYdWRELMqfpwHfAd5S5Px2P+P8FVB3WhpJkyStX1kH3gXMrznscuBwJdOBpyNiaZvraWYtKHirvknlcVpeRns0N5aPS7pD0plVb/cUepRXwL8CV0n6mKRTgFOBwnmI2tqrHhFPAZeNccirgUvTiCMmAudFxFWSjsnnnwpcSRqKtJA0HKmlpEpm1gHFZkd6vHa4YTVJ1wKvGWXXZ0mjab5Iatt+EfgP4KOl61lHRFyd485s4HFgp4h4pOj53c6rfj+jNIVzwKysB3BsN+tlZsW1o/MHICLG6gtZdT3pO8BP8se2PMqT9C/AB4E9gR2A6yR9MiKuKHJ+Pw5HMrN+F2q8tKBmCOIhrHqkdzkwU9LakrYijfe+uYlLbAxMi4jfRsS3gX2BE4qe7GnlzKy8zncO/bukHfOVFgF/DxARd0q6CLgLWA4c20yPekScUPP5QWCfouc7cJpZaerwDPAR8ZEx9p0CnNJMuZL+MyJOkPRjRgn/EXFQkXIcOM2snDY94+yR7+WfXx7zqAYcOM2svAENnBExV9IEYFZENJ3Rz4FzyBxw2992pNxOZc98d4eyZ85Zcl5Hyp0658MdKXfgDGjgBIiIEUmvL/OmUC0HTjMrbYBv1SvuJ6XPuBxYOdlFRHylyMkOnGZW3uAHzvvysgawft5W+Ldy4DSzcga7c6jiroj4QfUGSR8oerIHwJtZeSsKLP3tMwW3jcotTjMrRQxui1PS/qS5MLaQ9PWqXRuQBtQX4sBpZuUNaOAEHgbmAAcBc6u2PwucWLQQB04zK2eAn3FGxO3A7ZLOi4hlzZbT9WeckiZLuljS3ZIWSHprzX5J+nqeGv8OSTt3u45m1sDgP+PcMsehuyTdX1mKntyLzqGvAVdFxF+TpphbULO/OsvlLNK8fGbWRzqdc6gLvkuKLcuBdwDnAN8venJXA6ekDUnz350BEBEvRcSfag5zlkuzfhcFlv62bkT8DFBEPJjTAx1Y9ORutzi3Av4IfFfSbZJOzyk0qjnLpVk/izQ7UqOlz70oaQ3gXkkfl3QI8IqiJ3c7cE4Edga+FRE7kV51+nQzBUXEaRExNSKmTtywNvaaWUcNfovzeGA94DhgF+AjwBFFT+52r/piYHFE3JQ/X8zqgdNZLs363AA8wxxTRNySV/9ME3nNup1z6BFJf5C0bUTcA8wgzeRc7XJSdrsLgF1xlkuz/jOggTNP6lFXP09k/A/AuZLWIs1QcpSzXJoNjgHpNa/nraQ+lPOBm0gvQpXW9cAZEfOA2pShznJpNkgGN3C+hpRb6FDgw8AVwPkRcWeZQjzJh5mVNqjjOCNiJCKuiogjgOmkO9vrJH28TDl+5dLMyuv/4UZ1SVqbNGbzUGBL4OvApWXKcOA0s3L6uEXZiKRzgO1JfSknR8T8BqeMyoHTzMob0MAJHEYaP348cJy0sm9IpC6WDYoU4mecZlZap59xSvqApDslrZA0tWbfZ/IkQPdI2rdq+35520JJo75YExFrRMT6edmgalm/aNAEtzitoE5lz/xxh7JnHtCh7Jm3PXxBR8oF2OmWmR0ru9268ErlfOC9wLdfdl1pO2Am8GZgc+BaSW/Mu79J6jFfDNwi6fKIqB0n3hYOnGZWThdeqYyIBQBVt9IVBwMXRMSLwAOSFgLT8r6FEXF/Pu+CfGxHAqdv1c2svGLvqm9SmYgnL7PacOV6kwAVmhyoXdziNLNSROFb9ccjovZll1XlSNeSBqTX+mxEXNZc7brDgdPMSlO0fq8eEXs3cdpYkwB1bXIg36qbWTlFbtM79wz0cmCmpLUlbUXKFHEzcAuwjaSt8jwYM/OxHeEWp5mV1ukB8Hli4W8AmwJXSJoXEftGxJ2SLiJ1+iwHjo2IkXzOx4GrgQnAmWXfPy/DgdPMSuv0cKSIuJQ6r0FGxCnAKaNsv5L0RlDHdTvn0LaS5lUtz0g6oeYYZ7k063eDPwN8S7o9kfE9wI4AkiaQHt7W/qtSneVyV1Imul27WE0zG8sAv6veLr3sHJoB3BcRD9Zsd5ZLsz5WGY404MnaWtLLwDmTNAtzLWe5NOt3EY2XIdaTwJmHCxwE/KDZMpzl0qx3BnUi43bpVa/6/sCtEfHoKPuc5dKsn42Dzp9GenWrfiij36ZDGrR6eO5dn46zXJr1HY00XoZZ11uckiaRpn76+6ptznJpNkCG/Va8kV5kuXwO2Lhmm7Ncmg2KYOg7fxrxm0NmVtqwDzdqxIHTzEoRvlV34DSzcsbBOM1GHDjNrDS3OM3MSvIzTjOzMgJYMb6bnA6c1lPv6VDa4cs6lHZ43807k3YYOpN6eNp6T7a9TGDcvznkwGlmpcktTjOzctw5ZGZWhif5cOA0s3LSAPjxHTkdOM2sNI04cJqZFedb9e7PxynpREl3Spov6XxJ69TsX1vShTnL5U2Stux2Hc1sLAXSZgz5rXy30wNvARwHTI2I7UmJ42fWHHY08FREvAH4KvClbtbRzBrTimi4DLNezAA/EVhX0kRgPeDhmv0HA2fn9YuBGZLUxfqZ2Vii81kuJX0g35mukDS1avuWkv4iaV5eTq3at4uk3+W71a93Mm50NXBGxBLgy8BDwFJSWoxrag5bmeUyIpYDT1Mz8bGZ9Vjnb9XnA+8Fbhhl330RsWNejqna/i3g74Bt8rJfq5Wop9u36huRWpRbAZsDkyQd1mRZTg9s1itRYGml+IgFEXFP0eMlbQZsEBE35iwS5wB/01ot6uv2rfrewAMR8ceIWAb8EHhbzTErs1zm2/kNgSdqC3J6YLPe0YoVDRdgk0rjJi+z2nT5rSTdJul6SXvkbVsAi6uOWZy3dUS3hyM9BEyXtB7wF2AGMKfmmMuBI4DfAu8Hfp7/BTGzfhBAsWeYj0fE1Ho7JV0LvGaUXZ+NiMvqnLYUeF1EPCFpF+BHkt5cqDZt1NXAGRE3SboYuBVYDtwGnCbpC8CciLgcOAP4nqSFwJOs3utuZj0koi1vDkXE3k2c8yLwYl6fK+k+4I2kO9UpVYdOyds6ohdZLj8PfL5m8+eq9r8AfKCrlTKzclb0ZiZjSZsCT0bEiKStSZ1A90fEk5KekTQduAk4HPhGp+rRi+FIZjbIKrfqjZYWSDpE0mLgrcAVkq7Ou/YE7pA0jzRc8ZiIqEw6+jHgdGAhcB/w09ZqUZ9fuTSz0jo9yUdEXApcOsr2S4BL6pwzB9i+oxXLHDjNrLxx3l/rwGlm5UT07Blnv3DgNLPyxnfcdOA0s/I8kbHZEDq4Q9kzf/hgZ7JnQmcyaP4+VnvprnUBjIzvJqcDp5mVNPzzbTbiwGlm5TlwmpmV5MBpZlZCBIyM9LoWPeXAaWblucVpZlZCAEOeU6gRB04zK2+cvznUi/TAx+fUwHdKOmGU/cqJlhZKukPSzt2uo5mNxemBu51zaHtSMqVpwFuAd0t6Q81h+7Mq2dIsUgImM+sXQWpxNlqGWLdbnG8CboqI53MGy+tJmeyqHQycE8mNwOSciMnM+oVbnF01H9hD0sY579AB5MRsVVamB85GTbrkLJdmvRLplctGyxDrds6hBZK+BFwDPAfMA5oaEBYRpwGnAay3zebD/c+bWT8JiBjuwNhI1zuHIuKMiNglIvYEngJ+X3PIyvTAWUeTLplZE1ZE42WI9aJX/VX55+tIzzfPqznkcuDw3Ls+HXg6IpZ2uZpmNpZx/oyzF+M4L5G0MbAMODYi/iTpGICIOBW4kvTscyHwPHBUD+poZvX4lcuepAfeY5Rtp1atB3BsVytlZqXEkA83asRvDplZScN/K96IA6eZlRP4Vr3XFTCzwRJADHmveSMOnGZWTgSM83GcDpxmVtp4b3EqhuAhr6Q/Ag8WPHwT4PEOVGPQyu1k2S63f8p9fURs2s6LS7oq16GRxyNiv3Zeu18MReAsQ9KciJg63svtZNkudzDLteK6/uaQmdmgc+A0MytpPAbO01xux8t2uYNZrhU07p5xmpm1ajy2OM3MWuLAaWZW0lAGTklnSnpM0vw6+5vKpCnptZJ+IemunHi5kuQAAAReSURBVKXz+HaULWkdSTdLuj2Xe/Iox6wt6cJc7k2StixS53zuBEm3SfpJm8tdJOl3kuZJmjPK/ma/58mSLpZ0t6QFkt7aarmSts31rCzP1GZZbaG+J+Y/t/mSzpe0Ts3+pr5jOSNs/4qIoVuAPYGdgfl19h8A/BQQMJ2UQK5IuZsBO+f19Umz12/Xatn52Ffk9TWBm4DpNcd8DDg1r88ELizxfXyCNGH0T0bZ10q5i4BNxtjf7Pd8NvC3eX0tYHI7yq06fwLwCGlweKt/dlsADwDr5s8XAUe2+h0D25NydK1HesPvWuAN7fwevDS/DGWLMyJuAJ4c45CmMmlGxNKIuDWvPwssYPVEcqXLzsf+OX9cMy+1vXYHkwIKwMXADElqVGdJU4ADgdPrHNJUuQWV/i4kbUj6h+8MgIh4KSL+1Gq5NWYA90VE7dtmzZY7EVhX0kRSoHt4lHLLfsfOCNvHhjJwFlAok+ZY8u3WTqTWYctl59vpecBjwOyIqFtu/ov0NLBxgar+J/ApoN6sDM2WCym4XyNprqRZY5WdFfkutgL+CHw3P144XdKkNpRbbSZwfjvqGxFLgC8DDwFLSalerqlXbonvuG0ZYa39xmvgbImkVwCXACdExDPtKDMiRiJiR1JyummStm+1TEnvBh6LiLktV3B0u0fEzsD+wLGS9mxDmRNJj1m+FRE7kbKhfroN5QIgaS3gIOAHbSpvI1LLbytgc2CSpMNaLTciFgCVjLBX0UJGWGu/8Ro4m86kKWlNUtA8NyJ+2M6yAfJt6S+A2skRVpabbwk3BJ5oUNxuwEGSFgEXAO+U9P02lFup65L88zHgUmBavbKzIt/FYmBxVYv7YlIgbbXciv2BWyPi0VH2NVPu3sADEfHHiFgG/BB4W71yy3zH4YywfWu8Bs6mMmnm51JnAAsi4ivtKlvSppIm5/V1gX2Au0cp94i8/n7g5xEx5tsLEfGZiJgSEVuSbk9/HhG1raHS5eZ6TpK0fmUdeBfp9rK27FLfRUQ8AvxB0rZ50wzgrlbLrXIoo9+mN1vuQ8B0Sevl/z9mkJ5915bbzHfsjLD9qte9U51YSH8xlpIyaS4GjgaOAY7J+wV8E7gP+B0wtWC5u5Oe691BunWaR3r21FLZwA7Abbnc+cDn8vYvAAfl9XVIt5cLgZuBrUt+J3uRe9XbUS6wNXB7Xu4EPpu3t+N73hGYk7+PHwEbtancSaSW3oZV29pR7smkf+jmA98D1m7Td/xL0j8atwMz2lVfL60vfuXSzKyk8XqrbmbWNAdOM7OSHDjNzEpy4DQzK8mB08ysJAdOaytJJ0nyUA0bah6OZG2VJxWZEmnSCbOh5MBpZlaSb9WtrXyrbuOBA6eZWUkOnGZmJTlwmpmV5MBpZlaSA6eZWUkOnGZmJTlwmpmV5MBpZlaS3xwyMyvJLU4zs5IcOM3MSnLgNDMryYHTzKwkB04zs5IcOM3MSnLgNDMryYHTzKyk/w/SnR1iJUkl6wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "IOGoOa35zHCk" }, "source": [ "$\\mathbf{y}$ is the unknown vector which is contains the numerical approximations of the $y$. \n", "\\begin{equation}\n", "\\color{red}{\\mathbf{y}}=\\color{red}{\n", "\\left(\\begin{array}{c} y_1\\\\\n", "y_2\\\\\n", "y_3\\\\\n", ".\\\\\n", ".\\\\\n", "y_8\\\\\n", "y_9\n", "\\end{array}\\right).}\n", "\\end{equation}" ] }, { "cell_type": "code", "metadata": { "id": "imGFAqN4zHCk" }, "source": [ "y=np.zeros((N+1))\n", "# Boundary Condition\n", "y[0]=1\n", "y[N]=2" ], "execution_count": 4, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "PlD_mMG6zHCl" }, "source": [ "and the known right hand side is a known $9\\times 1$ vector with the boundary conditions\n", "\\begin{equation}\n", "\\mathbf{b}=\\left(\\begin{array}{c}-99+3x_1^2\\\\\n", "3x_2^2\\\\\n", "3x_3^2\\\\\n", ".\\\\\n", ".\\\\\n", "3x_8^2\\\\\n", "-209+3x_9^2 \\end{array}\\right)\n", "\\end{equation}\n" ] }, { "cell_type": "code", "metadata": { "id": "xvJAKYbKzHCm", "outputId": "26f32aaa-008e-4a3d-9331-9b064ee3a521", "colab": { "base_uri": "https://localhost:8080/" } }, "source": [ "b=np.zeros(N-1)\n", "for i in range (0,N-1):\n", " b[i]=3*h*(i+1)*h*(i+1)\n", "# Boundary Condition\n", "b[0]=-y[0]*(1/(h*h)-2*1*h/(2*h))+b[0]\n", "b[N-2]=-y[N]*(1/(h*h)+2*9*h/(2*h))+b[N-2]\n", "print('b=')\n", "print(b)" ], "execution_count": 5, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "b=\n", "[-9.8970e+01 1.2000e-01 2.7000e-01 4.8000e-01 7.5000e-01 1.0800e+00\n", " 1.4700e+00 1.9200e+00 -2.1557e+02]\n" ] } ] }, { "cell_type": "markdown", "metadata": { "id": "PqkEWxHmzHCm" }, "source": [ "## Solving the system\n", "To solve invert the matrix $A$ such that \n", "\\begin{equation}A^{-1}Ay=A^{-1}b\\end{equation}\n", "\\begin{equation}y=A^{-1}b\\end{equation}\n", "The plot below shows the graphical representation of $A^{-1}$." ] }, { "cell_type": "code", "metadata": { "id": "GmWj2LbJzHCm", "outputId": "3a61f449-c561-4acb-c000-debb9746cf5a", "colab": { "base_uri": "https://localhost:8080/", "height": 297 } }, "source": [ "invA=np.linalg.inv(A)\n", "\n", "plt.imshow(invA)\n", "plt.xlabel('i',fontsize=16)\n", "plt.ylabel('j',fontsize=16)\n", "plt.yticks(np.arange(N-1), np.arange(1,N-0.9,1))\n", "plt.xticks(np.arange(N-1), np.arange(1,N-0.9,1))\n", "clb=plt.colorbar()\n", "clb.set_label('Matrix value')\n", "plt.title(r'Matrix $A^{-1}$',fontsize=32)\n", "plt.tight_layout()\n", "plt.subplots_adjust()\n", "plt.show()\n", "\n", "\n", "y[1:N]=np.dot(invA,b)" ], "execution_count": 6, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVUAAAEYCAYAAADsymWcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3debxdVX338c+XjCRAEhKGSBgrD2ipZYgB61A0jFpBq7baqkGl6KNWQPtYqX0Mg1qqPGqpLYoQAQcccGBQwYBS2yrUgMyoAURICGMigwyZfs8fe51k53Du3evcu/c5N/d+36/Xet199l577XXvTX537Wn9FBGYmVk9tuh3B8zMRhMHVTOzGjmompnVyEHVrIck7SXpbyWdJ+lmSWslhaST+903q8f4fnfAbIz538Dx/e6ENccjVbPeugU4A/hr4HnAl/rbHaubR6pmPRQR55Q/S1rfr75YMzxSta6k638h6e5+98VsJHJQbSPp7lLgCEmPSprSxf4ntu0fkt7VZJ9Lxz45lRN6cbzNiaTT2n4nf9fvPtno5NP/atsArwcuyKz/tgb7UmVh+vpb4DN97MeIIml3oD2I7tuPvtjo55Hq4Frv8GYFSkkvBP4offS1spHj08Bk4OnSOgdVa4RHqoP7ETAf+FNJe0TEXRX1356+rgeuBl7RYN/6IiLU7z50Q9KhwNHp46nA3wPTgL0kTY6IpwfceWMb/wS8dgiHPykivjOE/Wwz5qA6uEUUQVXAMcBHBqooaTLwxvTxSmBF052zwUmaAPxL+ngPxYj1z4A/ofi3vw+wJKOp2cBeQ+jCtCHsY5s5n/4P7nrgprS8QNJgP6/XAdPT8hdzGlfhRZJOkbRY0r2SnpL0tKT7JF0h6QRJWw3Sxm6tmy+l1bt2uFkWkq7utF9rvaQJkt4m6XJJ90hanbbvVtpvwLv/kvZJ/Q9J6yQNOlKXNK90jDWSDsr5uXXhbymeBYVi1Pg0xXOiLVmXACLimIjQEMp5NX8/thlwUK22KH3dBThkkHqt666rgNxTvp8AP6UYAR8CzKG49jeJYnR0GMXo6k5JL+uu292RtEvqyyLgcGBnYEI3bUTELcAH0sctgC9JmjnA8bYGvlo6xsKIuGYIXe9I0vZsPLO4FrgwLZeD6n51Hc+sxaf/1b4MfAKYSBE4f9heQdKubLx++tWIeEbKuvS4Xfq6nOI//lLgUYpA8wfAEcD2qfxA0gsj4ra2NlYC/yctfzJ9XQV8vMPx7h2gH5OB71IEmQeBy4DfAFsBL2bjDbtKEfHvkg6juI75HIpR+1Edqv47xfcIxfXn03OPkel0Np5+fyA2Thx8c6mOb1ZZ/SLCpVSAuymCSAB7p3UXpc9PATM67HNyaZ8D0rrzSuveNcCxTgFeOEhfJlIE9FY7V1X0vVXv7ozvc7dS/VY5D5g63GMAM4Flpbrvbdv+ltK2h4Gdav4dvpDiZmEA32zbNqt07McB9fnfW+vfycn97IdLjb/TfndgpJUBguorS+ve3VZfpX1uLK2vDKpd9OniUlvPHaTecILqT3ICTO4xgIOBdaU/Ri9I658LPFZq5zU1//4EXJPafgbYo0Od+0vH37PH/772T/1rlYdSP5a1rZ/dy3651Fd8TTXPFcB9afntbdvmA7um5awbVENwYWn5pQ0d49RI/+vrEBFXA/+UPk4GLpQ0jeI66tZp/eci4rt1HTNZAByYlv81Oj8G1/XNqhptQ9G/VpmV1u/Utn5Sj/tlNfE11QwRsU7SBcCHgAMkvSAiWk8FtILsaorrr0MiaSfgBRTXIbdm09/NH5SWh/JoT5WngB830O7JFNeaXwQ8n+JJil3SttuA99d5MEnbsDGQPwJ8dICqN1P8MYQiqH6zzn4MJv2x2aye9bXuOKjm+yJFUIXihtWJkqaz8aHwSyPi4W4blfSXFDeaDsjcZXp1la4tjYh1dTcaEWsl/RVwA8VNo1ZAfRp4Y0Q8VfMhFwI7puVTIuJ3A9TzEwDWGJ/+Z4qIXwP/nT6+OT1Y/iaKU1vo8tQ/PaN6DvA18gMqpePV6dEG2gQgIu5m4+ix5eSIuLlD9SGTtDfFc6kAvwLOGqR6P0//bZTzSLU7iygeMZoFvJqNp/4rgMu7bOsdqQCsSW1fCtxO8VjTU63Ro6SXU7wy25TG5ilIo/l3t60+UtInI6LO436Gjc+8/h44Z5DH2rYsLc+WtH1EPFhjX2wMc1DtzjeAM4GpwGkU1wkBLhjC6fP7SstviIiLB6m7Ob/u+AU2nvYHxfXEPwVOAj5WxwEkHU3xwkLL/qnk2pcOzx+bDYVP/7sQEU+w8abG80ubuj31n8LG2azuqgioULyjvtmRdCzFtIlQPMb0KopROcDJdbyWKmkS8KlhNuNLAFYbB9XuLWr7/NOI+FWXbcwoLa/MqP+6zHbXpq/juutO/STtxcY5XQNYEBE/AP4xrRsPfDXdsR+OvwP2SMtfjcz38imeTGhxULXa+PS/SxHxn5I+CbSyAVSNMjtZxcZT4edJ2nKgO+GSjiH/P/2jFG8zzaiq2CRJEymerZ2aVn0qIlqn158EDqWY62B3itdV3zzE48yhuIwAxYP+H+5i99tLy34CwGrjkeoQRMQHI+K9qSwewv5PsnHKuanA59Np7CZSQP08+e/et0bMUyXN67ZfNTqdjYHqeuAfWhvSCwZvpXg9FeCvJb1liMf5JBsD92fTkwa5ykF1T0lbDljTrAseqfbP6cC30vJbgJdJ+j7F64qzKGao+kOKVz0/Tt4o7HsUc4UCXCrpKxSpVVo30ZZHw5MmSzocaOXI+j3wVxGxulwnIlZIejtwSVr1b5J+GhF3dnGcl7Jx/tpVdH/T69cUP5dxqbyAYlIbs+Hp93uyI63Q4d3/IbZzXqmdgSZU+QgbJ/7oVJ6gyA9/cGndeYMcc2uKEdhA7V1dqrtbp/UZ39eA7/5TzKZVfq/+2Iq2/q1U91pgfGYftgB+Udr3/UP8Hf261MY7+/1vz2V0FJ/+91FEnErxLv83KKb/W0Nx4+omipHsCyLiK1209zjFe+P/CPwstbV20J1qouKh0POAHdKqb0VbjvsOPsDGB/HnUTymluOdbLzOfDfw2eyObqp8CcA3q6wWiqhtDg0zszHPI1Uzsxo5qJqZ1chB1cysRg6qZmY1clA1M6uRg6qZWY1GxRtVE6ZNiUk7DHdejmfbQs08bjZOzUxfOq6h/k7YovakAACMVzPtTmio3Ulq5pHfSQ393m68ac3DEbFddc1qh798ajyysvrnet1Nz1wREUfUcczN1agIqpN22IY/+uyC2tvdasLq6kpDMG1S3VlEUrsTnm6k3e0mPt5Iu9tPfKyRdp8zYVUj7e42oetsOVn2GN/Mv7Md5qz4bV1tPbxyHddeMaey3oTZd86qrDTKjYqgamZNC9bVmqhh9HJQNbNKAazPnixtbHNQNbNKQbCm/oS7o5KDqpll8Ug1j4OqmVUKYJ2DahYHVTPL4pFqnp4+/C9pkaQHJd0ywHZJOlPSHZJuktRNmmEza0gAayIqi/X+jarzgMEeDD4S2DOV44CzetAnM6sQBOsyivU4qEbETxg8JfPRwAVRuAaYLml2b3pnZgMKWJdRbORdU90JuLf0eVlat6K9oqTjKEazTNy+/ldUzWyjQKxB/e7GZmGznVAlIs6OiLkRMXfCNGcXNmtSAOujutjIG6kuB3YufZ6T1plZn63zSDXLSBupXgK8NT0FcBDwaEQ869TfzHqruPu/RWWxHo9UJV1IkcN+lqRlwEJgAkBEfA74PvBK4A7gSeBtveyfmXVWPPzvkWqOngbViHhTxfYA3tOj7phZpkCsG3EntiPTSLumamYj1PrwSDWHg6qZVQrE6hjX725sFhxUzaxSMZ+qT/9zOKiaWRbfqMrjoGpmlSLEGp/+ZxkVQXWcgmmT6k96t+yy3WpvE+CerRtpltXTmskhtH5aM1lEt5rxZCPt7jqjmcR/B8y4p5F2D9v65kbarVPxSJVP/3OMiqBqZk0T6/xwfxYHVTOr5BtV+RxUzaySH6nK56BqZlnW+/Q/i4OqmVXyjap8DqpmVinwI1W5/KfHzCpFwLrYorIMh6RtJS2WtDR9nTFAvQWpzlJJC0rrD5B0c0oceqYkpfUnS1ou6YZUXjmsjlbodTbVnSX9WNJtkm6VdHyHOs6oajbiiPUZZZg+BFwVEXsCV6XPm/ZC2pZiytADgXnAwlLwPQv4GzYmDy0nGf10ROybyveH29HB9Hqkuhb4QEQ8HzgIeI+k57fVcUZVsxEmgNUxvrIM09HA+Wn5fOA1HeocDiyOiJURsQpYDByREoRuExHXpClELxhg/8b1Opvqioi4Pi0/DtxOkdivzBlVzUaYQKyP6kIxAf2SUjmui8PsUMr0cT+wQ4c6AyUH3Sktt69veW8681000GWFuvTtRpWk3YD9gGvbNmVlVC1nU528Q0PvfZrZBpl3/x+OiLkDbZR0JbBjh00fLn+IiJBUVyrBs4DTKAbcpwH/D3h7TW0/S1+CqqStgG8BJ0TEY0NpIyLOBs4GmLbXDs7jaNagIpvq8E9sI+KQgbZJekDS7IhYkc5OH+xQbTlFSqaWOcDVaf2ctvXL0zEfKB3jC8BlQ+1/jp7f/Zc0gSKgfiUivt2hijOqmo0wrUeqqsowXQK07uYvAC7uUOcK4DBJM9Jp/GHAFemywWOSDkp3/d/a2r/t8uFrgVuG29HB9Pruv4Bzgdsj4lMDVHNGVbMRaB2qLMN0OnCopKXAIekzkuZKOgcgIlZSnML/PJVT0zqAdwPnUCQOvRP4QVr/ifSo1U3Ay4ETh9vRwfT69P/FwFuAmyXdkNb9A7ALOKOq2UgVocZfU42IR4D5HdYvAY4tfV4ELBqg3j4d1r+l3p4OrtfZVP8LBv9z5oyqZiNPgN+oyuTXVM0sg+dTzeWgamaVirv/zlGVw0HVzLJ4lqo8DqpmVikQa31NNYuDqplVKmap8ul/jlERVMdvsY5Zk5+ovd07ZjbzotZulzaTRfSp2ZMbafeJHSc20u6Tz5nQSLu37rRlI+0+s66Z/y67THykkXbhrlpb8zXVPKMiqJpZszxJdT4HVTOr5Lv/+RxUzSxD829UjRYOqmZWKQLWOKhmcVA1syweqeZxUDWzSq2Z/62ag6qZZakhsd+Y0NOgKmky8BNgUjr2RRGxsK3OJIqkXQcAjwB/GRF397KfZrapANau9yNVOXp9keQZ4BUR8cfAvhRZEA9qq/MOYFVEPBf4NPDPPe6jmbXLSPrnywOFXmdTjYhovfo0IZX215bKaWovAuanjAFm1idBcfpfVaw/OarGpVn/H6TI3z1gNtWIWAs8Cszs0M5xrTS4T//u6aa7bTamFaf/W1QW60NQjYh1EbEvRUK/eZKelf4gs52zI2JuRMydPL2Zd97NbCOf/ufp25+WiPgd8GPgiLZNG7KpShoPTKO4YWVmfdJ6pMpBtVqvs6luJ2l6Wt4SOBT4ZVu1cpra1wM/SnmrzKyPfE01T6+fU50NnC9pHEVA/0ZEXCbpVGBJRFxCkcL6S5LuAFYCb+xxH82sTQS+Zpqp19lUbwL267D+I6Xlp4E39LJfZlbNp/d5/EaVmVXya6r5HFTNLItTVOdxUDWzShE+/c/loGpmWcJBNYuDqpllEOt89z/LqAiqE7ZYz+zJj9Xe7ppZa2tvE+DRPac00u7M797aSLtT9t61kXYf3XNqM+2undRIu7/dckYj7f5m5naNtFsn56jK5z89ZlYtiuuqVWU4JG0rabGkpelrx79ikhakOkslLSit/5ikeyU90VZ/kqSvS7pD0rWSdhteTwfnoGpmWXrwRtWHgKsiYk/gqvR5E5K2BRYCBwLzgIWl4HtpWteup9OJOqiaWaVI11SryjCVp/08H3hNhzqHU8xutzIiVgGLSfOHRMQ1EbGiot3GpxN1UDWzLJmn/7NaU3KmclwXh9ihFBTvB3boUGfD1KDJsrRuMFnTiZZJeomkt6Xl7STtXt39wqi4UWVmzct8pOrhiJg70EZJVwI7dtj04U2PFSGpLxMpSVoIzAX2Ar5IMZn+l4EX5+zvoGpmlSKo5ZGqiDhkoG2SHpA0OyJWSJpNMZF9u+XAwaXPc4CrKw7bmk50WeZ0oq+lmKPk+tTn+yRtXXGMDXz6b2ZZmr77z6bTfi4ALu5Q5wrgMEkz0g2qw9K63HZzphNdnbYHgKSunv3rS1BNKVV+IemyDtt6+viDmeWJUGUZptOBQyUtBQ5Jn5E0V9I5RR9iJXAa8PNUTk3rkPQJScuAKZKWSTo5tXsuMDNNJ/p+OjxV0OYbkj4PTJf0N8CVwBdyv4l+nf4fD9wObNNh24bHHyS9keLxh7/sZefMbFNBLUFz8GNEPALM77B+CXBs6fMiYFGHeh8EPthhfVfTiUbEGZIOBR6juK76kYhYnLt/z4OqpDnAq4CPUfzVaHc0cHJavgj4rCR59n+zPhpjE6qkIJodSMv6cfr/GYq/JusH2N51NtUnVz3TVF/NrCUyyigg6XFJj6XytKR1krLfg+/pSFXSnwEPRsR1kg4eTlsRcTZwNsCOf7jtKPl1mo1cY2WWqojYcKc/vSRwNHBQ7v69Hqm+GDhK0t3A14BXSPpyWx1nUzUbYQJYv16VZbSJwncp3uTK0uscVScBJwGkkerfRcSb26q1Hn/4Gc6majYyBDBGRqqS/rz0cQuKFwGezt1/RDz872yqZiPfGBravLq0vBa4m+ISQJa+BdWIuJr0JoSzqZqNdCJG4el9JxHxtuHsPyJGqma2GRjlI1VJ/8og32VEvC+nHQdVM6sWY+Lu/5I6GnFQNbM8o3ykGhHnV9eq5qBqZnlG/0gVKOZPBf4eeD4wubU+Il6Rs79nqTKzPGPkjSrgKxRzk+wOnEJx9//nuTuPipHqRK1ll0n1vx+wzfZPVFcagsf2aCYr55Yv2auRdid9P/vfU1emr9+nkXbXTsme+rIrK7fbspF2Vzw9rZF2azWGnlMFZkbEuZKOj4j/AP5D0tgKqmbWvBhoto7RZ036ukLSq4D7gG1zd3ZQNbM8Y2ek+lFJ04APAP9KMUXpibk7O6iaWZb+ZIzqi2sj4lGKGfJe3u3OvlFlZtVCsD6jjA7/LemHkt6RUrZ0xUHVzPKMkbv/EfG/gH8E/hC4TtJlktonfhqQg6qZ5RkjQRUgIv4nIt4PzKOY2Cn7xQAHVTPLM0aCqqRtJC2Q9APgp8AKiuCapR85qu4GHgfWAWsjYm7bdgH/ArwSeBI4JiKu73U/zawkQKPnmmmVG4HvUmRq/Vm3O/fr7v/LI+LhAbYdCeyZyoHAWemrmfXTKBmJZthjOBPjj8TT/6OBC1Iag2socm/P7nenzGxsGG6mkQGDasogOC8tr0+fBysPSrpY0h5VfQZ+KOk6Scd12L4hm2qyLK1r79+GbKpPrFpd/Z2a2bBovSqLDX76fypFQGstV0XvbYDXUmQ4PWSQei+JiOWStgcWS/plRPwkt8Mt5Wyqu+yzzdg5MTHrh1F0I2ooJE2MiKzR24BBNSJOKS2fnHngq4ELB6sTEcvT1wclfYfirlo5qG7IpprMSevMrJ/GSFBNceyYiLg7fZ4HfAH445z9676m+l/AXw+0UdJUSVu3loHDgFvaql0CvFWFg4BHI2JFzf00sy4pqsso8U/A5ZLeLeljwOeA7LxVtd79j4hVwMWDVNkB+E7x1BTjga9GxOWS3pX2/xzwfYrHqe6geKRqWEm4zKwmY2SWqoi4IsWkxcDDwH4RcX/u/j19pCoi7qLDEDoF09ZyAO/pZb/MbHCjbCQ6KEn/F/gL4GXAC4CrJX0gIr6Xs/9IfKTKzEaiUHUZBknbSlosaWn62nEyk/S209JUFpTWf0zSvZKeaKt/jKSHJN2QyrEVXZkJzIuIn0XE54HDgRNyvw8HVTPLovXVZZg+BFwVEXsCV6XPm/ZB2hZYSPFC0DxgYSn4XsrAr5N+PSL2TeWcwToRESdExFOlz7+NiENzvwnPp2pmeZo//T8aODgtnw9cTZGAr+xwYHFErASQtBg4ArgwvSxEumfTNUmfiYgTJF1Kh+82Io7KacdB1cyq5V9TnSVpSenz2emZ8hw7lJ70uZ/ixna7rJeDOnidpJcBvwZOjIh7O9T5Uvp6RmZ/O3JQNbM8eaf3D7dPklQm6Upgxw6bPlz+EBEh1XZr7FKKkewzkt5JMQp+VrrpiLhO0jjguIgY8NHQKqMiqE7UWnaeUH821b1mPVh7mwBLdm8mK+dDz0yurjQEs8Zlz3rWlan/tbSRdqfstGcj7T722LhG2l21upl/D3WrI8RFxIBvW0p6QNLsiFiR5vvo9B9wORsvEUDxctDVFccsB4dzgE8MUnedpF27eYOq3agIqmY2KlwCLABOT187PfN+BfDx0s2pw4CTBmu0FajTx6OA2yv6cRdFSpVLgN+3VkbEpyq/A3z338xyNT9J9enAoZKWUswfcjqApLmSzgFIN6hOA36eyqmlm1afkLQMmCJpmaSTU7vvk3SrpBuB9wHHVPTjTuAyivi4dSpb5X4THqmaWbWo5ZGpwQ9RnKbP77B+CXBs6fMiYFGHeh8EPthh/UlUjGbb3BYR3yyvkPSG3J09UjWzPGMknQqdA3B2UPZI1cwqidH/mqqkIynmHdlJ0pmlTdsAa3PbcVA1s2o9OP0fAe4DllDczLqutP5x4MTcRhxUzSzPKB+pRsSNwI2SvhoRa4baTs+vqUqaLukiSb+UdLukF7Vtl6QzJd0h6SZJ+/e6j2bWwdi5prpbilG3SbqrVXJ37seNqn8BLo+IvSmmAWx/ZqycTfU4imyqZtZnY2iS6i9SxJ21wMuBC4Av5+7c06AqaRrFHIXnAkTE6oj4XVs1Z1M1G2mC4jXVqjI6bBkRVwFKM1SdDLwqd+dej1R3Bx4CvijpF5LOSWlVyrrOpvrYyuwbc2Y2RGNopPqMpC2ApZLeK+m1dPHwf6+D6nhgf+CsiNiP4hWwZ82ZmCMizo6IuRExd5ttfb/NrHFj55rq8cAUirevDgDeQvHabJZeR6NlwLKIuDZ9vohnB1VnUzUbgcbAI1UARMTP0+ITDCFHXq9zVN2f0h3sFRG/ongl7ba2apcA75X0NYrZvZ1N1azfRtdItKM0gcqARvIk1X8LfEXSRIrZYN7mbKpmI5tSGeVeRHE/50LgWob4Lfc8qEbEDUD7JLbOpmo2wo2B0/8dgUOBNwF/BXyPYnLrW7tpxBOqmFmeUX6jKiLWRcTlEbEAOIjibPlqSe/tph3fNjezPJt50MwhaRLFM6lvAnYDzgS+000bDqpmVm10PYfakaQLgH0o7uucEhG3DKUdB1UzyzIGrqm+meLZ+eMpsgW01ovids82OY04qJpZnlE+Uo2IWu4xjYqgOlFr2WX8qtrbPWDaPbW3CfDELpMaaffX4zulSR++Fds0k6V1mx33bqTdGbc/1Ui745+Y0ki7j69p5udbt9F++l+XURFUzaxhrQlVrJKDqplVGgvpVOrioGpmeRxUsziomlkWhaNqDgdVM6s2NhL/1cJB1czyeKCaxUHVzLL4RlWeXueo2kvSDaXymKQT2uo4m6rZSJNO/6uK9X6S6l8B+wJIGkcxo3/7ZAXlbKoHUmQ1PLCH3TSzTjxSzdLPqf/mA3dGxG/b1jubqtkI03pOdYwk/huWfgbVN1LMsN2u62yqq1b6vMOsaVoflWVY7UvbSlosaWn6OmOAegtSnaWSFqR1UyR9T9IvJd0q6fRS/UmSvp4uKV4rabdhdbRCX4JqSqVyFPDNobZRzqY6Y1vPtW3WqJwJqoc/Uv0QcFVE7AlcRYdMy5K2BRZSXBKcBywsBd8zImJvYD/gxZKOTOvfAayKiOcCnwb+edg9HUS/otGRwPUR8UCHbc6majYC9eBG1dHA+Wn5fOA1HeocDiyOiJURsQpYDBwREU9GxI8BImI1cD1F7Ghv9yJgvkrz+tWtX0H1TXQ+9Ycim+pb01MAB+FsqmYjQ95IdVbrslwqx3VxhB1K/9fvBzpNu1Z5eVDSdODVFKPdTfaJiLXAo8DMLvrVlZ4/pyppKkVyrXeW1jmbqtlIFuReM304ItoTe24g6UqKBHvtPrzJ4SJC6v7Wl6TxFAO2MyPirm73r0M/sqn+nra/EimYtpadTdVsBKrj7n5EHDJg+9IDkmZHxIr0xM+DHaotBw4ufZ4DXF36fDawNCI+07bPzsCyFHSnAY8M7Tuo5js8Zpan+RtVlwAL0vIC4OIOda4ADpM0I92gOiytQ9JHKQLmCW37lNt9PfCjNHhrhIOqmVVSVD9ONdxHqoDTgUMlLQUOSZ+RNFfSOQARsRI4Dfh5KqdGxEpJcyguITwfuD69sXlsavdcYKakO4D30+Gpgjr53X8zy9L0w/0R8QjFS0Ht65cAx5Y+LwIWtdVZRvGOQqd2nwbeUGtnB+GgamZ5/MZUFgdVM8vi11DzjIqgOkHwnPFra2937pRmnsiYtP2aRtqdveVjjbT7q5nbN9Lu/btmpVHv2qrnNZOddEpDT0uvWTeumYbrFMA6R9UcoyKomlnzPFLN46BqZnmcoyqLg6qZVXOOqmwOqmZWqZhP1SPVHA6qZpbHI9UsDqpmVi1/QpUxz0HVzDKEb1Rl6vm7/5JOTOkObpF0oaTJbdt7mvrAzPI4R1WeXqeo3gl4HzA3IvYBxlHkqirraeoDM8sUUV2sL7NUjQe2TPMaTgHua9ve09QHZpYhQOuisliPg2pELAfOAO4BVlCkSvlhW7Ws1AflbKqPPOLbkmaNa34+1VGh16f/MyhGorsDzwGmSnrzUNoqZ1OdOdPTwpo1TRGVxXp/+n8I8JuIeCgi1gDfBv6krc6GbKq9SH1gZhlaE6pUFet5UL0HOEjSlHSddD5we1udnqY+MLNqonqU6pFqoafPqUbEtZIuosjJvRb4BXC2pFOBJRFxCUXqgy+l1AcrefbTAWbWDw6aWfqRTXUhsLBt9UdK23ua+sDMMng+1Wx+o8rMsvj0Po+DqpnlcVDN4qBqZhn8xlQuB1Uzq+ZrqtkcVM0si6+p5hkVQXUCW7D9uKm1t7v/xMdrbxdp37oAAAh7SURBVBPgOeNua6Td/be8u5F27992WiPt3rfzjEbaXb53M+3e+1Qz7S57fHoj7dbOQTXLqAiqZtawCFjnOTZyOKiaWR6PVLM4qJpZHgfVLA6qZlYtAOeoyuI588wsQ8D6ddVlGCRtK2mxpKXpa8c7g5IWpDpLJS1I66ZI+p6kX6Z0TaeX6h8j6SFJN6Ry7LA6WsFB1cyqtUaqVWV4PgRcFRF7Alelz5uQtC3F3CEHAvOAhaXge0ZE7A3sB7xY0pGlXb8eEfumcs5wOzoYB1Uzy9N8jqpyKqXzgdd0qHM4sDgiVkbEKmAxcEREPBkRPy66GaspZsKbM9wODUU/sqkenzKp3irphA7bJenMlE31Jkn797qPZtYuYP366gKzWmmOUjmui4PsEBEr0vL9wA4d6mxIt5QsS+s2kDQdeDXFaLfldSmeXCRp5y761LWe3qiStA/wNxTD9tXA5ZIui4g7StWOBPZM5UDgrPTVzPolaAXNKg9HxNyBNkq6Etixw6YPb3K4iJC6T3qdsoVcCJwZEXel1ZcCF0bEM5LeSTEKfkW3befq9d3/5wHXRsSTAJL+A/hz4BOlOkcDF6TZ/q+RNF3S7NJfMDPrhxoeqYqIQwbaJumB1v91SbOBBztUWw4cXPo8B7i69PlsYGlEfKZ0zHI6pnPYNN7Urten/7cAL5U0U9IU4JWkfFQllcN72DSb6kOPDO+uo5lVSW9UVZXhKadSWgBc3KHOFcBhkmakG1SHpXVI+ihFTrtNLiumAN1yFM9O4VSrXqdTuV3SPwM/BH4P3AAMKSJGxNkUf5WY+8eT/QCdWZMCIhp/TfV04BuS3gH8FvgLAElzgXdFxLERsVLSacDP0z6npnVzKC4h/BK4vkiBx2fTnf73STqKIoXTSuCYJr+JfqRTOZciDxWSPk4xEi3bkE01mZPWmVk/NfzwfzpNn99h/RLg2NLnRcCitjrLAA3Q7knASbV2dhD9uPu/ffq6C8X11K+2VbkEeGt6CuAg4FFfTzUbAZp/pGpU6Mdrqt+SNBNYA7wnIn4n6V0AEfE54PsU11rvAJ4E3taHPppZWQSs872LHP04/X9ph3WfKy0H8J6edsrMKkXeI1VjnidUMbMMPr3P5aBqZtUCn/5nclA1s0oBhKf+y+KgambVIqD551RHBQdVM8vikWoexSi4+CzpIYo3MHLMAh5uoBtu1+2OtHZ3jYjt6jiopMvTsas8HBFH1HHMzdWoCKrdkLRksFl03K7bHYvtWn08SbWZWY0cVM3MajQWg+rZbtftul1rypi7pmpm1qSxOFI1M2uMg6qZWY1GZVCVtEjSg5JuGWD7kDK2StpZ0o8l3ZaywR5fR9uSJkv6H0k3pnZP6VBnkqSvp3avlbRbTp/TvuMk/ULSZXW1K+luSTdLukHSkg7bh/oznp4yXv5S0u2SXjTcdiXtlfrZKo+pLZPvMPp7Yvqd3SLpQkmT27Z3/fOVMw5v3iJi1BXgZcD+wC0DbH8l8AOKmcIPokhGmNPubGD/tLw18Gvg+cNtO9XdKi1PAK4FDmqr827gc2n5jcDXu/h5vJ9iMvDLOmwbUrvA3cCsQbYP9Wd8PnBsWp4ITK+j3dL+4yjSH+9aw+9tJ+A3wJbp8zeAY4bz8wX2ocjlNoXijccrgefW+TNwabaMypFqRPyEIhfNQDZkbI2Ia4Dp2jQ52EDtroiI69Py4xQJxNqTEnbddqr7RPo4IZX2O4hHUwQcgIuA+ZI6po8oU5G751UUWSQ7GVK7Gbr+OUiaRvEH8VyAiFgdEb8bbrtt5gN3RkT7G3hDbXc8sKWK1MhTgPs6tNvNz3dDxuGIWAu0Mg7X0VfrgVEZVDNkZWwdTDqN249iVDnsttMp+g0UaXkXR8SA7ab/bI8CMzO6+hngg8BAs2EMtd0AfijpOknHDdZukvNz2B14CPhiulxxjqSpNbRb9kaKvPDD7m9ELAfOAO4BVlCk/vnhQO1m/nxryzhs/TFWg+qwSNoK+BZwQkQ8VkebEbEuIvalSHQ4T9I+w21T0p8BD0bEdcPu4LO9JCL2B44E3iPpZTW0OZ7iss1ZEbEfRcbdD9XQLgCSJlKkKP5mTe3NoBg17g48B5gq6c3DaTMibgdaGYcvZxgZh60/xmpQHXLGVkkTKALqVyLi23W2DZBOd38MtE9KsaHddKo5DXikorkXA0dJuhv4GvAKSV+uod3WKI2IeBD4DjBvoHaTnJ/DMmBZaZR+EUWQHW67LUcC10fEAx22DaXdQ4DfRMRDEbEG+DbwJwO1m/vzjYhzI+KAiHgZsIri2v1w+2o9MlaD6pAytqZrYecCt0fEp+pqW9J2kqan5S2BQynyl7e3uyAtvx74UUQM+uZGRJwUEXMiYjeK094fRUT7SKrrdiVNlbR1axk4jOK0tb3drn4OEXE/cK+kvdKq+cBtw2235E10PvUfarv3AAdJmpL+bcynuM7e3m63P19nHN6c9ftOWROF4j/OCoqMrcuAdwDvAt6Vtgv4N+BO4GZgbma7L6G4lngTxWnZDRTXvIbVNvAC4Bep3VuAj6T1pwJHpeXJFKetdwD/A+zR5c/kYNLd/+G2C+wB3JjKrcCH0/o6fsb7AkvSz+K7wIya2p1KMUKcVlpXR7unUPwBvAX4EjCphp/vf1L8MbkRmF9XX116U/yaqplZjcbq6b+ZWSMcVM3MauSgamZWIwdVM7MaOaiamdXIQdVqI+lkSX6cxMY0P1JltUmTt8yJYpIPszHJQdXMrEY+/bfa+PTfzEHVzKxWDqpmZjVyUDUzq5GDqplZjRxUzcxq5KBqZlYjB1Uzsxo5qJqZ1chvVJmZ1cgjVTOzGjmompnVyEHVzKxGDqpmZjVyUDUzq5GDqplZjRxUzcxq5KBqZlaj/w9r8DTFYiX/KgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "tNaP_8KHzHCn" }, "source": [ "## Result \n", "The plot below shows the approximate solution of the Boundary Value Problem (blue v)." ] }, { "cell_type": "code", "metadata": { "id": "cjjHauFqzHCn", "outputId": "2341c73d-216d-4f6d-a8a3-023f9319a4db", "colab": { "base_uri": "https://localhost:8080/", "height": 279 } }, "source": [ "fig = plt.figure(figsize=(8,4))\n", "\n", "plt.plot(x,y,'v',label='Finite Difference')\n", "plt.xlabel('x')\n", "plt.ylabel('y')\n", "plt.legend(loc='best')\n", "plt.show()" ], "execution_count": 7, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAAEGCAYAAAB1pazcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAXeklEQVR4nO3df3BV5Z3H8c+Xy9VYJYIk6bAGDLumKkKIIaj4Y5sacSnTkdGFtSytytSFcWvD7HYQdZ1lR12ttnVkp7Ysax3WLQXxx1TW0q2rVWGxdBsQRUOXslYxLTsJ4I+CUgG/+8e9iQn5dUNy7jnPzfs1c8fcc57c8/Xh5n7uc85zzjF3FwAACM+wuAsAAADHhxAHACBQhDgAAIEixAEACBQhDgBAoIbHXUB/lZSUeEVFRdxlAACQN1u2bNnr7qXHLg8uxCsqKtTY2Bh3GQAA5I2ZvdXdcnanAwAQKEIcAIBAEeIAAAQquGPi3Tl8+LCam5t16NChuEvBICkqKlJ5ebnS6XTcpQBAYhVEiDc3N2vEiBGqqKiQmcVdDgbI3bVv3z41Nzdr/PjxcZcDAIlVECF+6NAhAryAmJlGjx6t1tbWuEsBhoSZyzaqac/7XZZPGFOs9YsujaEi5KpgjokT4IWFf08gf2rGjVQ61flvLp0y1ZwxKqaKkKuCGIkDQJIlfaTbUF+px7Y0S/rk1tQpMzXUnxlfUchJwYzE45ZKpVRdXd3+ePPNN3XRRRf1+Xs33HCDmpqaJEl33313v7dbV1ens846S1VVVTr77LN100036d13321f37GGxYsX69xzz9XixYvV2tqqCy64QOedd542btzY7+0CyF3SR7plxUWaM6W8vcZ0yjS7dqzKRhTFXNknZi7bqIpbftzlMXPZ0P78GnIj8ai+EZ900knatm1bp2UvvfRSn7/30EMPtf98991367bbbuv3tletWqXa2lp99NFHuvXWWzVr1iy9+OKLXWpYsWKF9u/fr1QqpTVr1mjSpEmdtt+Xo0ePKpVK9bs+YKgLYaTbscak1SZlvgj9uuX3Onz0kz5M0hehuAy5kXg+vxGfcsopkqQXXnhBdXV1mj17ts4++2zNmzdP7pk3Yl1dnRobG3XLLbfoww8/VHV1tebNmydJ+sEPfqDzzz9f1dXVWrhwoY4ePdrr9k444QTdd9992r17t1555ZVONVx55ZU6cOCApkyZonvvvVc333yznnrqKVVXV+vDDz/UM888o2nTpqmmpkZz5szRgQMHJGUuc7tkyRLV1NToscce67Xd0qVLVVNTo0mTJulXv/qVJOnAgQOaP3++Jk2apKqqKj3xxBOS1OPrAIUohJFuW41mSlxtUuZLxrBj5sok8ctGvg25EI/qjdAWwNXV1brqqqu6rH/55Zf1wAMPqKmpSW+88YY2bdrUaf03vvGN9tH8qlWrtGPHDj366KPatGmTtm3bplQqpVWrVvVZRyqV0uTJk9tDtM26devaX3/JkiW64447dM0112jbtm06ePCg7rrrLj377LPaunWramtrdf/997f/7ujRo7V161ZdfvnlvbYrKSnR1q1bdeONN+pb3/qWJOnOO+/Uqaeequ3bt+vVV1/VZZddpr179/b6OkAh6vjZk9Twaaiv1NSK0xJZWwhfhOIQ2e50Mxsr6RFJn1ZmH9IKd192TBuTtEzSTEkfSLre3bdGVZP0yRvh0ca3dfioD9obobvd6R2df/75Ki8vl6T2Y+aXXHJJj+2fe+45bdmyRVOnTpWU+ZJQVlaWUy1to/xcbd68WU1NTbr44oslSR999JGmTZvWvv6aa67Jqd3VV18tSZoyZYqefPJJSdKzzz6rNWvWtLcZNWqUnn766V5fByhEbZ89q/57d2LDp6y4SGsXJvdvMcm7/OOavBjlMfEjkr7u7lvNbISkLWb2n+7e1KHN5yVVZh8XSPpe9r+RiuONcOKJJ7b/nEqldOTIkV7bu7uuu+463XPPPf3aztGjR7V9+3adc845Of+Ou2v69OlavXp1t+tPPvnknNq1/T/29f/X1+sAhaqhvlI7Ww4kKnxCkuQvQnEds49sd7q772kbVbv77yXtkHT6Mc1mSXrEMzZLGmlmY6KqqU1Sj/2k02kdPnxYklRfX6/HH39cLS0tkqT9+/frrbe6vRNdu8OHD+vWW2/V2LFjVVVVlfN2L7zwQm3atEm7du2SJB08eFA7d+487nYdTZ8+XQ8++GD783feeee4XgfoSUizlttGukn5zAlRUnf5x3XMPi/HxM2sQtJ5kn5xzKrTJb3d4Xmzuga9zGyBmTWaWeNgXcUriW+EBQsWqKqqSvPmzdOECRN011136YorrlBVVZWmT5+uPXv2dPt78+bNU1VVlSZOnKiDBw/qqaee6td2S0tLtXLlSs2dO1dVVVWaNm1al2Pq/WnX0e2336533nlHEydO1OTJk/X8888f1+sAPUn66VsYXEn9IhTXMXvr7/HTfm/A7BRJL0r6R3d/8ph1T0v6hrv/V/b5c5KWuHtjT69XW1vrjY2dV+/YsaNfu48RBv5dkYuW9w/p0vue1x+OfNy+rGj4MG1Y8rnEfdCjsHV8Lw72e9DMtrh77bHLIx2Jm1la0hOSVh0b4Fm/lTS2w/Py7DIAyAmzlpEUcRyqjSzEszPPvy9ph7v3dP7QOknXWsaFkt5z9+73GQNAD0I4fQtDQ74P1UY5O/1iSV+WtN3M2s69uk3SOEly9+WS1itzetkuZU4xm3+8G3N3bppRQKI+zIPCkuRZyxha8n2aXmQhnj3O3WuqeuaT+qsD3VZRUZH27dun0aNHE+QFoO1+4kVFfBAjd5y+haGoIK6dXl5erubmZu4/XUCKioraL44D5CLpFyoBolAQIZ5OpzV+/Pi4ywAAIK+G3LXTAQAoFIQ4AACBIsQBAAgUIQ4AQKAIcQAAAlUQs9MBRCOueyQDyA0jcQA94g5hQLIR4gB6FNc9kgHkhhAH0CPuEAYkGyEOoFfcIQxILkIcQK/iuEcygNwwOx1An7hDGJBMhDiAPnGHMCCZ2J0OAECgCHEAAAJFiAMAEChCHACAQBHiAAAEihAHACBQhDgAAIEixAEACBQhDgBAoAhxAAACRYgDABAoQhwAgEAR4gAABIq7mAExmrlso5r2vN9l+YQxxVq/6NIYKgIQkshG4mb2sJm1mNlrPaw/1cz+3cxeMbPXzWx+VLUASVUzbqTSKeu0LJ0y1ZwxKqaKAIQkyt3pKyXN6GX9VyU1uftkSXWSvm1mJ0RYD5A4DfWVGmadQzxlpob6M2OqCEBIIgtxd98gaX9vTSSNMDOTdEq27ZGo6gGSqKy4SHOmlLePxtMp0+zasSobURRzZQBCEOfEtu9IOkfS7yRtl7TI3T+OsR4gFh1H44zCAfRHnCH+Z5K2SfojSdWSvmNmxd01NLMFZtZoZo2tra35rBGIXNto3EyMwgH0S5whPl/Sk56xS9JvJJ3dXUN3X+Hute5eW1pamtcigXxoqK/U1IrTGIUD6Jc4Q3y3pHpJMrNPSzpL0hsx1gPEpqy4SGsXTmMUDqBfIjtP3MxWKzPrvMTMmiUtlZSWJHdfLulOSSvNbLskk7TE3fdGVQ8AAIUmshB397l9rP+dpCui2j4AAIWOy64CABAoQhwAgEAR4gAABIoQBwAgUIQ4AACBIsQBAAgUIQ4AQKAIcQAAAkWIAwAQKEIcAIBAEeIAAASKEAcAIFCEOAAAgSLEAQAIFCEOAECgCHEAAAJFiAMAEChCHACAQBHiAAAEanjcBQBRmblso5r2vN9l+YQxxVq/6NIYKgKAwcVIHAWrZtxIpVPWaVk6Zao5Y1RMFQHA4CLEUbAa6is1zDqHeMpMDfVnxlQRAAwuQhwFq6y4SHOmlLePxtMp0+zasSobURRzZQAwOAhxFLSOo3FG4QAKDSGOgtY2GjcTo3AABYfZ6Sh4DfWV2tlygFE4gIJDiKPglRUXae3CaXGXAQCDjt3pAAAEihAHACBQkYW4mT1sZi1m9lovberMbJuZvW5mL0ZVCwAAhSjKkfhKSTN6WmlmIyV9V9KV7n6upDkR1gIAQMGJLMTdfYOk/b00+UtJT7r77mz7lqhqAQCgEMV5TPwzkkaZ2QtmtsXMro2xFgAAghPnKWbDJU2RVC/pJEk/N7PN7r7z2IZmtkDSAkkaN25cXosEACCp4hyJN0v6qbsfdPe9kjZImtxdQ3df4e617l5bWlqa1yIBAEiqOEP8KUmXmNlwM/uUpAsk7YixHgAAghLZ7nQzWy2pTlKJmTVLWiopLUnuvtzdd5jZf0h6VdLHkh5y9x5PRwMAAJ1FFuLuPjeHNt+U9M2oagAAoJBxxTYAAAJFiAMAEChCHACAQBHiAAAEihAHACBQhDgAAIEixAEACBQhDgBAoAhxAAAC1WeIm9nXzGxUPooBAAC5y2Uk/mlJvzSztWY2w8ws6qIAAEDf+gxxd79dUqWk70u6XtKvzexuM/uTiGsDAAC9yOmYuLu7pP/LPo5IGiXpcTO7L8LaAABAL/q8i5mZLZJ0raS9kh6StNjdD5vZMEm/lnRztCUCAIDu5HIr0tMkXe3ub3Vc6O4fm9kXoikLAAD0pc8Qd/elvazbMbjlAACAXHGeOAAAgcpldzrQrZnLNqppz/tdlk8YU6z1iy6NoSIAGFoYieO41YwbqXSq82UD0ilTzRlcGwgA8oEQx3FrqK/UsGOu/ZMyU0P9mTFVBABDCyGO41ZWXKQ5U8rbR+PplGl27ViVjSiKuTIAGBoIcQxIx9E4o3AAyC9CHAPSNho3E6NwAMgzZqdjwBrqK7Wz5QCjcADIM0IcA1ZWXKS1C6fFXQYADDnsTgcAIFCEOAAAgSLEAQAIFCEOAECgCHEAAAIVWYib2cNm1mJmr/XRbqqZHTGz2VHVAgBAIYpyJL5S0ozeGphZStK9kp6JsA4AAApSZCHu7hsk7e+j2dckPSGpJao6AAAoVLEdEzez0yVdJel7ObRdYGaNZtbY2toafXEAAAQgzoltD0ha4u4f99XQ3Ve4e62715aWluahNAAAki/Oy67WSlpjmTtglUiaaWZH3P1HMdYEAEAwYgtxdx/f9rOZrZT0NAEOAEDuIgtxM1stqU5SiZk1S1oqKS1J7r48qu0CADBURBbi7j63H22vj6oOAAAKFVdsAwAgUIQ4AACBIsQBAAgUIQ4AQKAIcQAAAkWIAwAQKEIcAIBAEeIAAASKEAcAIFCEOAAAgSLEAQAIFCEOAECgCHEAAAJFiAMAEChCHACAQBHiAAAEihAHACBQhDgAAIEixAEACBQhDgBAoAhxAAACNTzuAtC9mcs2qmnP+12WTxhTrPWLLo2hIgBA0jAST6iacSOVTlmnZemUqeaMUTFVBABIGkI8oRrqKzXMOod4ykwN9WfGVBEAIGkI8YQqKy7SnCnl7aPxdMo0u3asykYUxVwZACApCPEE6zgaZxQOADgWIZ5gbaNxMzEKBwB0wez0hGuor9TOlgOMwgEAXRDiCVdWXKS1C6fFXQYAIIEi251uZg+bWYuZvdbD+nlm9qqZbTezl8xsclS1AABQiKI8Jr5S0oxe1v9G0mfdfZKkOyWtiLAWAAAKTmS70919g5lV9LL+pQ5PN0sqj6oWAAAKUVJmp39F0k96WmlmC8ys0cwaW1tb81gWAADJFXuIm9nnlAnxJT21cfcV7l7r7rWlpaX5Kw4AgASLdXa6mVVJekjS5919X5y1AAAQmthG4mY2TtKTkr7s7jvjqgMAgFBFNhI3s9WS6iSVmFmzpKWS0pLk7ssl/b2k0ZK+a5lLix5x99qo6gEAoNBEOTt9bh/rb5B0Q1TbBwCg0MU+sQ0AABwfQhwAgEAR4gAABIoQBwAgUIQ4AACBIsQBAAgUIQ4AQKAIcQAAAkWIAwAQKEIcAIBAEeIAAASKEAcAIFCEOAAAgSLEAQAIFCEOAECgCHEAAAJFiAMAEChCHACAQBHiAAAEihAHACBQhDgAAIEixAEACBQhDgBAoAhxAAACRYgDABAoQhwAgEAR4gAABIoQBwAgUIQ4AACBGh7VC5vZw5K+IKnF3Sd2s94kLZM0U9IHkq53961R1dPRzGUb1bTn/S7LJ4wp1vpFl+ajBAAABizKkfhKSTN6Wf95SZXZxwJJ34uwlk5qxo1UOmWdlqVTppozRuWrBAAABiyyEHf3DZL299JklqRHPGOzpJFmNiaqejpqqK/UMOsc4ikzNdSfmY/NAwAwKOI8Jn66pLc7PG/OLuvCzBaYWaOZNba2tg54w2XFRZozpbx9NJ5OmWbXjlXZiKIBvzYAAPkSxMQ2d1/h7rXuXltaWjoor9lxNM4oHAAQojhD/LeSxnZ4Xp5dlhdto3EzMQoHAAQpzhBfJ+lay7hQ0nvuviefBTTUV2pqxWmMwgEAQYryFLPVkuoklZhZs6SlktKS5O7LJa1X5vSyXcqcYjY/qlp6UlZcpLULp+V7swAADIrIQtzd5/ax3iV9NartAwBQ6IKY2AYAALoixAEACBQhDgBAoAhxAAACZZn5ZeEws1ZJbw3iS5ZI2juIrzdU0Y8DRx8OHH04cPThwEXRh2e4e5ernQUX4oPNzBrdvTbuOkJHPw4cfThw9OHA0YcDl88+ZHc6AACBIsQBAAgUIS6tiLuAAkE/Dhx9OHD04cDRhwOXtz4c8sfEAQAIFSNxAAACRYgDABCoIRPiZjbDzP7HzHaZ2S3drD/RzB7Nrv+FmVXkv8pky6EP/9bMmszsVTN7zszOiKPOpOurHzu0+3MzczPjdJ9j5NKHZvYX2ffj62b2w3zXmHQ5/D2PM7Pnzezl7N/0zDjqTCoze9jMWszstR7Wm5n9U7Z/XzWzmkgKcfeCf0hKSfpfSX8s6QRJr0iacEybv5a0PPvzFyU9GnfdSXrk2Iefk/Sp7M830ofH14/ZdiMkbZC0WVJt3HUn6ZHje7FS0suSRmWfl8Vdd5IeOfbhCkk3Zn+eIOnNuOtO0kPSn0qqkfRaD+tnSvqJJJN0oaRfRFHHUBmJny9pl7u/4e4fSVojadYxbWZJ+tfsz49Lqjczy2ONSddnH7r78+7+QfbpZknlea4xBLm8FyXpTkn3SjqUz+ICkUsf/pWkB939HUly95Y815h0ufShSyrO/nyqpN/lsb7Ec/cNkvb30mSWpEc8Y7OkkWY2ZrDrGCohfrqktzs8b84u67aNux+R9J6k0XmpLgy59GFHX1HmWyg667Mfs7vdxrr7j/NZWEByeS9+RtJnzGyTmW02sxl5qy4MufThP0j6kpk1S1ov6Wv5Ka1g9Pcz87gMH+wXBMzsS5JqJX027lpCY2bDJN0v6fqYSwndcGV2qdcps0dog5lNcvd3Y60qLHMlrXT3b5vZNEn/ZmYT3f3juAvDJ4bKSPy3ksZ2eF6eXdZtGzMbrszuo315qS4MufShzOxySX8n6Up3/0OeagtJX/04QtJESS+Y2ZvKHEtbx+S2TnJ5LzZLWufuh939N5J2KhPqyMilD78iaa0kufvPJRUpc2MP5Canz8yBGioh/ktJlWY23sxOUGbi2rpj2qyTdF3259mSfubZ2QmQlEMfmtl5kv5ZmQDnGGT3eu1Hd3/P3UvcvcLdK5SZW3CluzfGU24i5fL3/CNlRuEysxJldq+/kc8iEy6XPtwtqV6SzOwcZUK8Na9Vhm2dpGuzs9QvlPSeu+8Z7I0Mid3p7n7EzG6S9FNlZmU+7O6vm9kdkhrdfZ2k7yuzu2iXMpMVvhhfxcmTYx9+U9Ipkh7Lzgnc7e5XxlZ0AuXYj+hFjn34U0lXmFmTpKOSFrs7e9aycuzDr0v6FzP7G2UmuV3PwOYTZrZamS+KJdl5A0slpSXJ3ZcrM49gpqRdkj6QND+SOvg3AQAgTENldzoAAAWHEAcAIFCEOAAAgSLEAQAIFCEOAECgCHEAAAJFiAMAEChCHECvzGxq9n7IRWZ2cvb+3BPjrgsAF3sBkAMzu0uZy26eJKnZ3e+JuSQAIsQB5CB7fe1fKnN/84vc/WjMJQEQu9MB5Ga0MtfFH6HMiBxAAjASB9AnM1snaY2k8ZLGuPtNMZcEQEPkLmYAjp+ZXSvpsLv/0MxSkl4ys8vc/Wdx1wYMdYzEAQAIFMfEAQAIFCEOAECgCHEAAAJFiAMAEChCHACAQBHiAAAEihAHACBQ/w8ObP34F3+AugAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] } ] }